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Preface 

In recent years, proposals for lunar and Mars missions have proliferated. A common feature 
of many of these proposals is a large aerobrake that would engender significant radiative heating 
during the atmospheric pass. These large acrobrakes are intended to perforin their maneuvers at 
high altitude; therefore, the flow may be in thermochemical nonequilibrium for all or at least a 
significant part of the trajectory. Existing equilibrium radiative heating codes are not well suited 
to this problem. Existing nonequilibrium radiative heating codes have been considered either too 
costly (NEQAIR, ref. 1) or too approximate (Viscous Shock Layer, RAD/EQUIL, ref. 2) for use in 
coupled solutions of this type of flow field. 

The objective of this work is to describe a new r method for efficiently predicting thermal radiation 
effects on vehicles flying through the rarefied upper layers of the atmosphere where thermochemical 
nonequilibrium flow occurs. 

The Langley Optimized RAdiative Nonequilibrium (LOR AN) code w^as developed for this 
purpose. The LOR, AN code includes atomic bound-bound, bound-free, and free-free transitions 
as w r ell as molecular band radiation. Radiative cooling effects can be treated if the LORAN code 
is coupled with a flow-field code. The Langley Aerotherrnodynamic Upwind Relaxation Algorithm 
(LAURA) code was used for that purpose in the development of LORAN. To achieve an efficient 
computation, spectral approximations and optimization are introduced and controlled by a few 
parameters in the code. To further reduce computation time, a radiation subgrid is selected. 

This user’s manual presents a review of the theory on which the LORAN code is based and is 
followed by a practical discussion of the implementation in a FORTR AN code. Descriptions of the 
input and output files and the run-time options are also provided. The manual should benefit both 
the occasional user of the LOR AN code and those who wash to use it as a starting point for further 
modeling enhancements. 

Part of the information presented in this report was included in a thesis entitled “Nonequilibrium 
Radiative Heating Prediction Method for Aeroassist Flowfields with Coupling to Flowfield Solvers" 
submitted by Lin C. Hartung in partial fulfillment of the requirements for the Degree of Doctor 
of Philosophy in Aerospace Engineering. North Carolina State University, Raleigh, North Carolina, 
April 1991. 
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THEORY 


1. Radiation Theory 

Predicting the emission and absorption of thermal radiation in a ga s involves computing the 
transitions between various energy levels in the atoms and molecules. The resulting complete 
spectral variation of emission and absorption is contained in the absorption coefficient k v and the 
emission coefficient f u . The absorption coefficient is obtained from the product of the population 
density N n of energy level n and its radiation absorption cross section o vn , summed over all the 
energy levels available. Thus, 

= (!) 

I 

where l is the absorbing levels. The spontaneous emission coefficient is found from a similar 
expression in which the amount of energy emitted hi/ u i also appears, and the absorption cross 
section is replaced by the transition probability for emission A u [ in which u and l denote the 
upper and lower energy levels, respectively. Thus, 

jr = E E hy ui A ui (2) 

U l 

where u is the emitting levels. These expressions and the remainder of the development in this 
section are written for a single species to simplify the notation. For a gas mixture such as air, 
properties are computed for each species and are added to obtain the mixture properties (i.e. , 

Kv — £, Kv S . 

In equilibrium gases, the electronic energy-level populations are determined as a function of a 
uniquely defined equilibrium temperature according to a Boltzmann distribution. The rotational 
and vibrational energy levels of each electronic state of a molecule are populated according to 
the equilibrium partition function. Equilibrium relations exist between species concentrations. 
All these equilibrium equations allow radiative properties to be referenced to temperature and to 
the concentrations of a few major species, thereby allowing the development of simple radiation 
step models. Also, under equilibrium conditions, the absorption and emission coefficients are 
related according to Kirchhoff’s law, so that equation (2) is not needed. 

In nonequilibrium gases, these simple relations for the energy-state populations and species 
concentrations no longer apply. The problem now is to determine the nonequilibrium populations 
of each energy level. The rotational and vibrational states in a molecule may be assumed 
populated according to distinct rotational and vibrational temperatures T r and T v , respectively, 
if the excitation/de-excitation mechanisms among the states of these energy modes are fast 
enough. If these mechanisms are not effectively instantaneous with respect to the t ime scales 
of the flow, further separation of energy modes is required. Generally, T r and T v adequately 
describe the respective energy modes in most hypersonic flows. The electronic states, however, 
being more widely separated, equilibrate more slowly, so that a single electronic temperature 
may be a bad assumption, even within a single species. In this situation the electronic energy- 
level populations in the absorption and emission coefficients must be calculated individually 
at every spatial and frequency point according to the specific local conditions. The emission 
coefficient cannot be calculated from the Kirchhoff relation between absorption and emission. 
The radiative cross sections, however, depend only on the configuration of an individual atom 
or molecule and are the same for equilibrium and nonequilibrium conditions. The necessary 
expressions for the nonequilibrium absorption and emission coefficients of the various radiative 
transitions are developed below without employing equilibrium relations. 




Most radiation work has historically been done in the Gaussian system of units. This system 
introduces the particular oddity of measuring the electron charge in statcoulombs. As this 
particular unit of measure is not commonly used in engineering, some of what follows may be 
confusing. To clarify the situation, relevant information on the Gaussian system is provided in 
appendix A. 

1.1. Atomic Mechanisms 

Atomic radiation involves transitions of electrons in bound or free states of an atom. 
Computing the radiation from atomic mechanisms requires knowledge of the energy levels of 
each atomic species. Many electronic energy levels occur in a typical atom and these must be 
modeled by a manageable set of inputs. To reduce bookkeeping and complexity, a common 
practice is to group neighboring energy levels into one lumped level along with an appropriate 
overall degeneracy. This practice has been adopted here. The grouped levels from the Park 
quasi-steady-state (QSS) method (ref. 1) are used for nitrogen and oxygen atoms (22 levels for 
atomic nitrogen, 19 levels for atomic oxygen), yielding values for the principal quantum number 
n and the energy level E n for each level. 

1.1.1. Bound-Bound Transitions 

Bound-bound or atomic line transitions are those between bound electronic levels of atomic 
species. The method used in the present study for computing bound-bound transitions is adapted 
from that developed by Nicolet (ref. 2) for the RAD/EQUIL code. This model resolves an atomic 
line by distributing an odd number of spectral points starting from the line center as shown in 
figure 1. Following the Nicolet development, with modifications required for nonequilibrium, the 
absorption coefficient in a line is given by 

7re 2 

— fi u Nib{v, N e , Tj>, T e , . . .) (3) 
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where //. u is the oscillator strength of the transition, TV/ is the population density of the lower 
state of the transition, and b is the line shape function that may depend on many variables. The 
constants e, m, and c are the electron charge, electron mass, and speed of light, respectively. 

The line shape is a normalized function that can be described by one of several theoretical 
models depending on the predominant broadening mechanism. The Lorentz shape describes 
lines in which Stark broadening by electron impact dominates, as is often the case for heavy 
atomic species. This line shape is given as 


7 S / 7r 

0' - VCL) 2 + (7 5 ) 2 


(4) 


where the dependence of 6 on the gas conditions has been omitted. Here, 7^ is the Stark half- 
width at half- maximum (HWHM) and ucl the frequency of the line center. The HWHM is 
the half-width of the line at half its height, a common descriptor for line shapes. 

Stark HWHM’s have been calculated by various authors, most notably Griem. (See ref. 3.) 
Page et al. (ref. 4) rioted that these line widths vary with a power n of the electron kinetic 
temperature T e as 

7 S (T e ) * (^)V(7o) = (^)”7o (Wc (5) 

thus allowing use of a single input reference for Stark broadening 7^, which is the HWHM per 
electron at the reference temperature Tq (10 000 K here). The assumption of a linear dependence 
on the electron number density N e is implicit in equation (5) and follows reference 5. The values 
of 7 q and n for the important lines treated in the LORAN code are included in the input files. 

The mechanism of resonance broadening yields Lorentz line shapes as well. This mechanism 
may exceed Stark broadening for situations with low' ionization. Resonance broadening occurs 
when radiating and perturbing particles are of the same type and uncharged (so Stark broadening 
does not occur). When the perturbers have the same natural frequency as the radiators, a 
resonance effect produces strong coupling between them and results iri a wader frequency spread 
(i.e., resonance broadening). This effect is important only when sufficient particles are in the 
appropriate initial state, generally when the initial state has a low-lying energy level. The 
RAD/EQUIL code, in fact, included this contribution only for a few lines with low-lying lower 
states. For completeness in nonequilibriurn heating, the LORAN code computes it for all lines. 
Hunt and Sibulkin (ref. 5) give the approximate HWHM as 1 
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where N a is the number of perturbing atoms per unit volume (i.e., the number of atoms in the 
lower state of the line transition) and gi and g u are the degeneracies of the low r er and upper 
states, respectively. The effective Lorentz line width is the sum of 7^ and 7^. (See ref. 6.) 

The second line shape is the Doppler profile, which describes line broadening due to thermal 
motion of the atoms. This Gaussian line shape is given as 


K = 


1 /hL2 

7Vv exp 


iy_ - t^CL ) 2 In 2 
“ (7 D ) 2 


(7) 


1 The 1/2 power on gi/g u was omitted in Nicolet. (See ref. 2.) Also note that the RAD/EQUIL code contains an extra 
gi factor, which results from a different treatment of N a . 
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The Doppler HWHM is given by 


D = -Cl (8) 

C V m 8 

where k is Boltzmann’s constant and m s is the mass of the radiating species s. The temperature 
is the heavy-particle translational temperature T% that is due to the thermal motion of the atoms. 

The Lorentz and Doppler line shapes are different, so that these two line widths cannot 
simply be added. Convolution leads to a Voigt line shape combining the two types of effects 
but requiring significantly more computation. Faster approximations to the Voigt line shape 
have been developed in recent years. The LORAN code applies an approximation developed by 
Whiting (ref. 7) if this additional accuracy is desired. This approximation gives the Voigt profile 
within 5 percent everywhere. The Voigt line width is computed from 7 ^ and 7 ^ according to 
reference 8 . If more approximate results are acceptable, both Lorentz and Doppler line shapes 
are computed, and the largest of b L and b D at each frequency point is used. 

After the appropriate line shape is calculated at each spatial grid point, the absorption 
coefficient in a line can be computed from equation (3) and from input data for the oscillator 
strength and line width and the predicted nonequilibrium energy-level populations. 

The total radiative emission from an atomic line into all directions (ref. 9) is given by 


J-Ul — N u^ulAul 


( 9 ) 


where A u [ is the Einstein coefficient for spontaneous emission of the transition from the upper 
atomic energy level u to the lower level l. Detailed balancing relates A u [ to the absorption 
oscillator strength f[ u (ref. 9) as follows: 


A u l = 


8n z e 
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(10) 


Combining these two equations leads to the spectral emission coefficient for bound-bound atomic 
transitions (now per steradian). The line shape b is the same as that for absorption. 
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When the electronic state populations are in equilibrium, the expressions for the bound-bound 
absorption and emission coefficients satisfy Kirchhoff’s law, the appropriate limiting form. 


LI. 2. Bound- Free Transitions 


Atomic bound-free or photoionization transitions involve a bound electron being detached 
from an atom by an incoming photon or the reverse process of radiative electron capture. A 
model for bound-free transitions can be developed by assuming a hydrogen-like atom. A good 
source of such a development for an equilibrium gas is Zeldovich and Raizer. (See ref. 9, ch. 5, 
sec. 5.) The radiative cross sections developed there also apply to nonequilibrium because the 
configuration of an atom in a particular excited state is unchanged whether the gas around it is 
in equilibrium or nonequilibrium. From that reference, therefore, the cross section for an atomic 
energy level n is 


1 32? i 4 mZ 4 e 10 Qf 
vn \/3 3^ch 6 n^g n 


( 12 ) 
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where Z is the charge number (Z = 1 for neutral species) and h is Planck's constant. Note that 
in this expression the electronic partition function of the ion has not been assumed equal to 1, 
as was done in reference 9. Further, the degeneracy has not been replaced in the denominator 
by using the hydrogenic model. For grouped or lumped electronic levels, the hydrogenie model 
does not correctly describe the degeneracy. 

Once the cross section is obtained, the spectral absorption coefficient can be computed for 
the bound-free process, given the nonequilibrium populations of the electronic energy levels of 
the atom. The relation is 

oc 

Kv = N n^vn ( 13 ) 

n* 

The lower limit n* indicates the photoionization threshold and is the lowest atomic energy level 
from which a photon of frequency v can detach an electron. The upper limit in practice has a 
finite value corresponding to ionization of the atom. 

Although quite accurate, equation (12) can be improved further. In particular, a correction 
from quantum mechanics (ref. 9, p. 266), sometimes called a Gaunt factor, may be applied to 
it. The approximate factor developed by Biberman (ref. 10) for low-energy photons has been 
applied here. Although only an average correction that was derived for thermal equilibrium 
conditions, the approximate factor nevertheless provides reasonable results for the flow-field 
conditions examined so far. More' detailed Gaunt factors may be available from other sources 
which apply better to nonequilibrium flows. However, these factors require substantially more 
work for a marginal improvement in the results because the order of magnitude' of the correction 
does not change. Therefore, inclusion of such detailed Gaunt factors has been deferred. 

Also significant is the error in this expression for low quantum states near the ground atomic 
state when' the hydrogen-like assumption is invalid. For many atoms, the cross section of the 
ground state is well-known from experimental data. This information is used for transitions 
involving the ground state instead of the value predicted by equation (12). 

Bound-free emission results from the capture of a free electron into an ion. The cross section 
for this process is given by Zeldovich and Raizcr (ref. 9) as 


128tt 1 ZV () 

3 \/3 uh"'’ h 1 v 2 n ,! v 


(M) 


where v is the initial speed of the captured electron. The number u r of such captures into the 
nth energy level for electrons with speeds between v and v + dv per unit volume per unit time 
is 

n c = N f N ( f ( v) dv va rn ( 1 5) 

The electron speeds follow the distribution function /(v). which is assumed to be a Maxwell 
distribution in equilibrium at the electron translational temperature T ( . (The validity of this 
assumption in a plasma with Debye effects may be a candidate for further investigation.) The 
energy emitted hu H [, where the subscript u now denotes the free state, is given bv the sum of 
the initial kinetic energy of the electron and the net energy required to reionize it from the final 
energy level /. 

2 

hv ul = ~ + (/ - E,) (16) 

where I is the ground-state ionization potential and Ej is the energy of the bound level and 
is measured from the ground atomic state. This equation is solved for v 2 and substituted to 
eliminate v in favor of the frequency v in equations (14) and (15). 
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After these pieces are combined, the emission coefficient for bound-free transitions is found 
from 

fv dv = ^ ^hv ul n c dv (17) 

where l is now the final levels. The complete expression is 


it 


128tt 4 Z 4 e w / m \3/2 n ^ l 

3v/3 m*<?h* N+Ne \2^k T e ) 


exp[—(his — I + E[)/kT e ] 


(18) 


The lower limit on the summation corresponds to capture of a zero-energy electron; n mm is 
determined for each frequency by setting v to zero in equation (16). The upper limit is 
determined by ionization of the atom. A quantum correction to this result is needed as discussed 
above. The Gaunt factor used for emission is identical to that for absorption. 

When the electronic state populations are in equilibrium with the free electrons, the 
nonequilibrium absorption and emission coefficients developed here for bound-free transitions 
satisfy Kirchhoff’s law, as required. 


1.1.3. Free- Free Transitions 


Atomic free-free transitions, also called bremsstrahlung radiation, result when a free electron 
passing near a heavy particle is decelerated and gives off energy in the form of radiation. Such 
radiation cannot occur in electron-electron collisions because momentum and energy cannot 
both be conserved. (See ref. 6, p. 137.) In fact, only transitions due to a nearby atomic ion are 
considered here. Free- free transitions caused by neutral atoms and molecules are assumed to be 
negligible because the cross sections are generally one or two orders of magnitude less. (See ref. 9, 
p. 256.) A model for free-free transitions can also be developed using semiclassical methods. 
Zeldovich and Raizcr (ref. 9, ch. 5, sec. 3) is again a good source. The result they obtained can 
be rederived without using equilibrium assumptions except for a Maxwell distribution of the 
energies of the free electrons at T e . In this case, the only change in the result is that T e replaces 
T. The equation for the free-free absorption coefficient can therefore be written as 

4/ 2t r \l/2 Z 2 e G KT XT 
K " _ 3 \3mkT e ) haw 3 + e 

where N+ is the concentration of the appropriate ion. This result could also be modified by a 
quantum correction, but the effect is minimal for the conditions in typical hypersonic flow fields. 

Free-free transitions are a special case in nonequilibrium gases because they involve only 
free electrons (the effect on the energy of the ion is negligible). Free electrons are assumed 
to equilibrate rapidly to some T e \ thus, Kirchhoff’s law and the development of the emission 
coefficient lead to the same result, 

,, 8 / 2tt y/2 z 2 e 6 

h= 3{3^rJ 1^ N + Nr (20) 

so that as T e — > T, the equilibrium values will be recovered. 

1.2. Molecular Mechanisms 

The treatment of molecular radiation requires consideration of vibrational and rotational 
transitions within the molecule as well as electronic transitions and combinations of the different 
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types of transition. Each such transition contributes a discrete line to the radiation spectrum 
whose line center frequency is determined by the energy differential of the transition. A line- 
by-line calculation of molecular radiation is possible, but it does not fit the requirements of 
the present study for a relatively rapid calculation method which can be coupled to a How-field 
code. As an alternative, therefore, the “smeared band” model of molecular radiation is used. 
In this model, the rotational line structure of the molecular radiation is smoothed to provide a 
continuous variation of the absorption coefficient. Because the rotational energy levels are so 
closely spaced, the lines resulting from rotational transitions are nearly overlapping except in 
rarefied conditions. The error introduced by this approximation should therefore be acceptable 
under most conditions of interest. However, if the LOR AN code is applied to a new gas regime, 
this assumption should be checked. 

The basic development of the smeared-band model can be found in several sources. (See 
refs. 9, 11, and 12.) The development given in appendix B is an enhancement of that found in 
Zeldovich and Raizer. (See ref. 9, ch. 5, sec. 3.) It incorporates higher order expressions for the 
molecular energy levels from Herzberg (ref. 13), while eliminating equilibrium assumptions. 

The final result for the smeared-band molecular absorption coefficient k is 
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The emission coefficient is given by 
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( 22 ) 


These expressions give the smeared absorption and emission coefficients for all rotational and 
vibrational levels in a given transition from lower electronic state A to upper electronic state B. 


1.3. Induced Emission 


The above expressions for the emission coefficient do not include all the radiation emitted 
by the medium. Stimulated or induced emission also occurs when photons are present. This 
emission is proportional to the radiative intensity l v and is. therefore, commonly included as a 
correction to the absorption coefficient. The induced emission (ref. 9. ch. 2. sec. 4) is given by 


7 __ r 
3\j Jv 


2/m 3 J 


( 23 ) 
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In equilibrium gases, I v in equation (23) becomes the Planck black-body function, and the usual 
result from Kirchhoff ’s law is obtained. In nonequilibrium gases, the intensity cannot be replaced 
with the Planck function. Instead, the induced emission coefficient of equation (23) must be 
substituted in the transport equation along with the coefficients for spontaneous emission and 
absorption. The transport equation then becomes 

+ «„(*)/„(*, fi) = jZ(8)+j' v (s) (24) 

where is the spontaneous emission coefficient and Q is the direction. (See sec. 6.) Substituting 
equation (23) and gathering terms containing I v then allows the definition of a corrected 
absorption coefficient for thermochemical nonequilibrium equal to the factor that multiplies 
I v . The result is 

c 2 

K v ~ K v — j u 2f iv ^ (25) 

The corrected absorption coefficient k! u may be negative in nonequilibrium. Negative values 
occur when the nonequilibrium populations are such that the second term on the right of 
equation (25) is larger than the first term. Physically, this negative k! v means that the upper 
energy level of a transition or transitions is overpopulated relative to the lower energy level. 
Such a population inversion can occur in a nonequilibrium boundary layer, for instance, when the 
higher energy levels are still populated according to the temperature and chemistry of the inviseid 
shock layer. This fact must be carefully considered when developing any solution algorithm for 
radiative transport in nonequilibrium problems. Because of the term 1 /j/ 3 in equation (25), 
negative values of k! u generally occur at the low-energy end of the spectrum. 

2. Electronic Excitation Theory 

The electronic states of atoms and molecules in a gas can become perturbed through many 
excitation mechanisms and reactions. In equilibrium gases, these effects balance, so that the 
states are populated according to a Boltzmann distribution. To obtain the nonequilibrium 
electronic-level population densities Nj for the ith electronic state of each atomic and molecular 
species in a gas, rate equations must be written for each level (of each species), including all 
possible populating and depopulating mechanisms. 

The populating mechanisms are as follows: 

1. Electron impact excitation or de-excitation from bound state j to bound state i: 

Nj + e~ — ► Nj + e~ Reaction rate: (KjjNjN e ) 

2. Heavy particle impact excitation or de-excitation from bound state j to bound state i : 

Nj -f- M — » Nj + M Reaction rate: (HjjNjNy) 

3. Radiative absorption in lower bound state j transitioning to higher bound state i : 

Nj 4* hv — ► N{ Reaction rate: (BjjU u Nj) 

4. Spontaneous radiative emission from higher bound state j decaying to lower bound state i : 
Nj — ► Ni + hv 
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Reaction rate: (AjjNj) 



5. Induced radiative emission from higher bound state j decaying to lower bound state i: 

Nj + hui — ► Ni + hvi + hu e Reaction rate: ( BjiU u Nj ) 

6. Collisional deionization from free state c to bound state i: 

7V+ + e~ + M — » Ni + M Reaction rate: f K cr N + N^) 

7. Radiative deionization from free state c to bound state i: 

N+ + e~ —* Ni + hu Reaction rate: ( A ci N + Np ) 

The corresponding depopulating mechanisms are as follows: 

1. Electron impact de-excitation or excitation from bound state i to bound state j: 

/Vj + er —> Nj + e~ Reaction rate: ( KijNiN e ) 

2. Heavy particle impact de-excitation or excitation from bound state i to bound state j: 

Ni + M -* Nj + M Reaction rate: (H^N^Nj,) 

3. Radiative absorption in lower bound state i into higher bound state j : 

Ni + hu — > Nj Reaction rate: {B l:j u u N r ) 

4. Spontaneous radiative emission from higher bound state i to lower bound state j: 

N t — ► Nj + fw Reaction rate: ( 

5. Induced radiative emission from higher bound state i to lower bound state j: 

Ni + hi', —> Nj + hu t + hu c Reaction rate: (BjjU^N,) 

6. Collisional ionization from bound state i to free state c: 

Ni + M -*N+ + e~ + M Reaction rate: ( K ir NiN c ) 

7. Radiative ionization from bound state i to free state c: 

Ni + hv — <■ N+ + e~ Reaction rate: ( Ai c u u Ni ) 

The rate of change of the population density Ni in a control volume can be written 

^ = N"' - N“ nt - V • (NiVi) (26) 

at 

where the first term on the right denotes the creation of through the first seven mechanisms 
above, the second term is the removal of N± due to the corresponding depopulation mechanisms, 
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and the last term represents the net flux of across the boundaries of the control volume. The 
velocity is the effective velocity of the particles in the 2 th electronic level of species s. The 
source and sink terms can be written explicitly in terms of the above mechanisms as 


Jmax jfmax i— 1 Jrnax 

N? = E K H N 3 N e + E HjiNjN h + BjiUvNj + J2 A ji N j 
3 = 1 3 = 1 3 = 1 j=i + 1 

Jmax 

T ^ ^ BjiUi/Nj T K^N+N^ + A^N-^Ne (27) 

j=i + 1 


Jmax Jmax Jmax i — 1 

^r it = E K a N i**e + E H a N i N h + E + E 

j=i j=i j=*+i j — i 

i— 1 

+ ^ B{jU v N{ -f Ki c N{N e + A{ c u v Ni (28) 

j=l 


To solve for the population densities, data on the various rates of excitation and de-excitation 
of each level are required. Radiative transition rates are reasonably well known. Collisional rates 
are more elusive. Sources for such data are varied, and the rates presented are often uncertain 
due to the experimental difficulties involved in the measurements. Park (ref. 1) has collected a 
comprehensive set of rates for nitrogen and oxygen species. He cites the uncertainty in the data 
set as a factor of 2 and asserts that the impact on the excited-state populations is negligible. 
This set has been adopted here; however, further verification of this accuracy would be desirable 
in the future. 


Equation (26) is not usually solved as is. Rather, the usual approximation is to assume a 
quasi-steady-state (QSS) distribution of excited states. The QSS model assumes that the rate 
of change in a population is small compared with the rate of transitions into and out of that 
level (ref. 14, ch. 3); that is, 

< Nf (29) 

< N° ut (30) 

If the divergence term V • (iVjVj) is also negligible, then equation (26) reduces to 


dNj 

dt 

dN x 

dt 


Nf x - N‘ mt « 0 


(31) 


or 



(32) 


and thus is transformed from a set of differential equations to a set of algebraic equations. This 
method is designed for use when the amount of nonequilibrium is “not too large,” such that the 
assumptions above are reasonable. In particular, the present application does not account for the 
history of the flow in areas of high gradients (boundary layer and shock), where y ■ (A^v*) is not 
necessarily negligible. Effectively, this approach decouples the excitation from the macroscopic 
flow-field solution. 
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The QSS assumption requires that the relaxation time be much less than the dwell time in 
the flow and is best satisfied for upper excited states in which the atoms and molecules transition 
rapidly back and forth to the free state. However, because QSS in effect predicts instantaneous 
accommodation to the local conditions, it will rise faster in the shock and fall faster in the 
boundary layer. The integrated effect on radiative flux to the wall, for example, may be quite 
accurate. 

A further difficulty arises in this equation because the radiative absorption and induced 
emission terms depend on the incident radiation in the form of radiation density u /y , which is 
not known a priori. A rigorous solution for the excited-state populations would require complete 
coupling of the radiative, excitation, and chemical processes in the flow, which is not vet feasible 
because of computational limitations and lack of rate data. An alternative would be an iterative 
process in which QSS populations without radiation are used to compute radiation, and the 
result is then fed back into equation (26) until the iteration converges. This process is formidable 
and computationally expensive, however. Instead, approximations are normally invoked. The 
simplest approximation is to assume that absorption and induced emission do not affect the 
QSS populations (optically thin assumption). This assumption is not accurate if ultraviolet lines 
arc considered, although induced emission can usually be neglected. To account for radiative' 
absorption, a radiation escape factor can be defined. (See refs. 14 16.) This approach is an 
approximate and somewhat empirical means of correcting the populations for absorption effects. 
It allows the limits on population levels to be determined by computing them with full absorption 
and no absorption. The escape factor is used here, based on an approximate flow dimension 
which the user must supply. The solution to the QSS equation proceeds according to the met hods 
developed by Park in the NEQAIR code. (See ref. 1.) 


3. Transport Theory 


The complete equation of radiative transport has been derived in many standard texts. 
Neglecting transients and assuming a nonscattering medium, it simplifies to 


rfJ„(.s,fi) 

ds 


+ = &{*) 


(33) 


However, neglecting scattering may not be valid for ultraviolet radiation in O 2 . for example; and 
scattering must be included to treat ablating flows. 

For some frequencies, n r u can be negative in part of a nonequilibrium flow (sec. 1.3). so it may 
pass t hrough zero somewhere in the domain. Therefore, the usual procedure of dividing by 
to transform this equation to an optical variable would introduce a singularity. For this work, 
then, equation (33) is left in terms of the physical coordinate. 

Equation (33) applies at each // and for each lb A formal solution can be written using 
an integrating factor. The radiative transport can then be computed numerically for a finite 
number of frequency points and directions. A method of selecting appropriate frequencies at 
which to compute the radiative properties Jf v and h J ly in a nonequilibrium flow is discussed in 
sections 5 and 6. Approximate methods for dealing with the directional variations are discussed 
next. 


3.1. Plane-Parallel Medium 

A common approximation for solving the radiative transport equation is to assume that 
the properties of the medium through which the radiation travels vary in only one direction. 
This situation is illustrated by the inset in figure 2 and may be a reasonable assumption in the 
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stagnation region of a flow field, where gradients are mainly in the normal direction. In this 
case s and H are replaced by 2 and /i (the cosine of the angle between Q and the 2 -axis). The 
derivative d/ds then becomes fid/dz. To accommodate the boundary conditions in this case, I u 
is split into forward (// > 0) and backward (// < 0) components. The formal solution for each 
can be written as 



(34) 


(z./i) = /„ ( 20 , n) exp 



r Is" 


+ 




(35) 


The boundary z — 0 is the vehicle surface, which is assumed to emit radiation according to 
Planck’s function at the imposed wall temperature T w . The boundary z = 20 I s the free-stream 
gas ahead of the vehicle, which is assumed not to emit radiation. The equations for the geometry 
along the normal to the wall (fi = ±1) are 


/+( 2 , 1 ) = I vp {T w )e~ T ‘'^ + ds' 

J 0 

/-(2, -1) = r o fu{ s')e^y ds > 


(36) 

(37) 
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where t v is the optical variable (introduced without division by 0 in equation (33)): 


Tv 



K»{ S ') ds' 


(38) 


To obtain the radiative flux in an absorbing, plane-parallel medium, the formal solutions 
above for the intensity must be integrated over all directions. If azimuthally symmetric radiation 
is assumed and if the exponential approximation is employed to replace the exponential integral 
functions (refs. 2 and 17), the result is 

f u (.s')e 2Tl ' {:) ~ 2Tl/{s,) ds' 

(39) 

Although this approach involves an approximation, it provides excellent results (within 10 per- 
cent for a typical blunt body) compared with the exact integration over the direction of the 
incoming intensity. 

The divergence of the radiative flux is obtained by differentiating equation (39) with respect to 
then integrating over the complete spectral range. The differentiation can be done analytically 
or numerically from the solution for q^(^). The latter approach has been implemented here and 
follows that of Nicolet. (See ref. 2.) The analytical approach requires a numerical differentiation 
of the integrand and would be similarly accurate. 

To obtain the total radiative flux, the spectral radiation is integrated over frequency: 


Qvr(z) = » 


lupiT^e- 2 ^ + 2 f Z jU^e 2 ^ 
J 0 


2t u(z) ds > _ 2 


L 


fV 2 

Qr(z) = / q v R {z)dv 
J V\ 


(40) 


The theoretical limits on the spectral integration are v i = 0 and 1/2 = oo. In practice, finite limits 
are required and are included as inputs to the computer code. The lower limit is often set around 
0.3 eV because a zero-energy photon causes a singularity in the free-free absorption coefficient. 
(See eq. (19).) The singularity arises because additional quantum mechanical considerations 
not discussed here must be included at very low energies. The radiative energy omitted by 
ignoring these low-energy transitions is minimal. Most gas radiation in shock layers in air 
occurs at energies below 16.5 eV, which is a typical upper limit. An alternate upper limit at 
6.2 eV includes only the visible and infrared portions of the spectrum. This range of energies is 
easier to measure experimentally and includes minimal self-absorption, but it is not sufficient if 
coupling effects are of interest. 

3.2. Multidimensional Medium 

When the tangent slab approximation is invalid, the directional variation of the radiation 
must be accounted for. Because an exact solution of the radiative transport equations in 
multidimensional media is not feasible, an approximate method is needed. One such method 
is described in reference 18. This method appears as an option in the LORAN code, but the 
interested reader is referred to that reference for more details because the transport method is 
still being developed. Overall, the required inputs and the resulting output are the same for 
either transport method. 

4. Flow-Field Coupling 

Radiation is inherently coupled to the state of the gas from which it emanates and results from 
radiative transitions in the internal energy states of the atoms, molecules, and ions in the gas. 
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The radiation emitted by such transitions travels through the gas at the speed of light until it 
is reabsorbed or exits the flow field. Radiation thus provides an additional energy transport 
mode. Radiation can also affect the chemical state of the gas through photodissociation, 
photoionization, and radiative recombination. 

Energy transport through the gas by radiation is modeled in the Navier-Stokes equations with 
a term that appears in the energy equation: the divergence of the radiative flux y * QR (denoted 
Qrad by Gnoffo, Gupta, and Shinn in ref. 19 ). This divergence term is computed in the LORAN 
radiation-transport model. When the tangent slab approximation is used, y • qp = dqp/c/77. The 
radiation effects on chemical kinetics can be included by adding photoreactions to the chemistry 
model of the gas and accounting for their effects on radiation properties. The radiation effects 
on chemical kinetics are not included in the current LAURA and LORAN codes. 

Energy transport due to radiation is usually neglected in the aerothermodynamic analysis of 
hypersonic flow fields. For flow fields in which this assumption is not valid, the effects of radiation 
can be determined by comparing solutions from calculations with and without radiative energy 
transport. In general, the shock layer becomes nonadiabatic when radiation transport is included. 
The enthalpy of a parcel of fluid decreases due to radiative cooling as it traverses the layer from 
the shock to the boundary layer edge. For equilibrium gas calculations, the temperature profiles 
decrease across the shock layer as a result. Lower average temperatures result in higher average 
densities and thus decreased shock standoff distance. The combination of lower temperatures 
and a thinner shock layer reduces the radiative heat load to the wall. In the boundary layer, 
however, a net absorption of radiation may occur, leading to locally increased temperature and 
enthalpy and to a possible increase in convective heating. 

In a thermochemical nonequilibrium gas calculation, multiple temperatures are defined to 
describe the partitioning of energy across the internal modes. The number of temperatures 
required to adequately describe the flow depends on the conditions. For example, rotational, 
vibrational, and electronic energy modes for the mixture or for each species may be used. 
The various states of each energy mode are assumed to be in equilibrium at the characteristic 
temperature of that mode. Each mode is described by its own energy conservation equation, 
which includes terms coupling it to the other energy modes. Usually, separate internal modes 
for each species are not defined because the energy modes are assumed strongly coupled between 
species. 

To correctly model the radiative coupling in the nonequilibrium case, only that portion of 
V • qR associated with transitions of each particular energy mode should be included in the ener- 
gy equation for that mode. Purely rotational radiative transitions occur at low infrared 
energies which are generally neglected in gas radiation — so all radiative transitions of interest in 
air are assumed to involve a transition between vibrational and/or electronic states. The LAURA 
code currently employs the two- temperature model. (See ref. 20.) This model is implemented 
with two energy equations: a total energy equation (for total temperature Tt) and a combined 
vibrational-electron-electronic energy equation (for the temperature Ty). Thus, because the 
complete v ■ qp contains no pure rotational transitions, it appears in both energy equations. 

PRACTICE 

5. Implementation 

The implementation in a practical computer code of the theory for radiation transitions and 
radiative transport developed in the section entitled “Theory” requires some additional work. 
The details of the process are described below. 

The present method is implemented in the LORAN code, which was developed using two 
existing codes as starting points. The RAD/EQUIL program (ref. 21) was followed in developing 
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some of the code structure. It also is the source of part of the algorithm used to model atomic 
line radiation. Significant changes have been made, however, to remove equilibrium assumptions. 
The QSS excitation portion of the NEQAIR program (ref. 1) has been minimally adapted for 
use in the LORAN code. 

The LORAN code is divided into four major sections: reading data and setting up properties 
which are constant for a calculation (including the optimized spectral array), calculating 
the nonequilibrium excitation, evaluating the radiative properties, and solving the radiative 
transport problem. The excitation calculation is from the well-known NEQAIR code; therefore, 
it will not be discussed in this section. 

5.1. Initial Setup 

5. LI. Data Input 

To evaluate the expressions for absorption and emission obtained in section 1 for each radiat- 
ing mechanism, input values are required. Data are also required to predict the rionequilibrium 
populations of the various bound energy levels; also, flow-field data are required to define the 
local gas conditions. Finally, program control inputs are desirable to maximize the flexibility of 
the method without requiring constant reprogramming and recompiling. These inputs should 
allow the user to turn individual radiating mechanisms or species on or off to perform detailed 
studies of the radiation spectrum. They should also provide for selection of options in the code 
(i.e., equilibrium vs. QSS populations, flux vs. intensity, etc). The inputs are obtained from var- 
ious sources, as discussed below. Most of the inputs are species properties and will not change 
unless significantly different, new data become available. 

5.1.2 . Flow- Field Data 

The information required about the gas conditions consists of the number densities of each 
radiating species considered and the temperatures for the several energy modes. The LAURA 
flow-field solver (ref. 22) is the source of the flow-field data used in developing this method. It 
is an 11-species model which provides number densities for all the nitrogen/oxygen radiating 
species of interest. The LAURA code currently incorporates a two- temperature model for 
thermal nonequilibrium. (See ref. 20.) The rotational and heavy-particle translational energy 
modes are assumed equilibrated, whereas the vibrational, electronic, and electron translational 
energy modes are assumed in equilibrium with each other but not necessarily with the other 
modes. Thus, the separate T v and T e carried in the development of the absorption and emission 
coefficients are equal, as are T r and T). The distinction between these temperatures was 
maintained in the LORAN code, however, so that each can be used if the separate temperatures 
become available from the LAURA code or from another CFD code in the future. Some recent 
work (ref. 23) suggests further that distinct vibrational temperatures may exist for each molecular 
species. Should a CFD code ever include such detail, minor modifications to the LORAN code 
will allow these individual values of T v to be used. 

An effort has been made to keep the input of the flow-field properties from the LAURA 
code generic so that, with minor changes to a single driver routine, the LORAN code can be 
interfaced with a different CFD code. In fact, a few cases have already been run with flow- 
field data obtained from the RAD/EQUIL code (ref. 2), a direct simulation Monte Carlo code 
(ref. 24), and a four-temperature CFD code. (See ref. 25.) Coupled solutions have been obtained 
with the LAURA and GIANTS codes. (See ref. 26.) 
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hv, eV 

Figure 3. Example of spectrum optimization for bound- free continuum. 

5.1.3. Maximum Rotational Quantum Number 

The radiative property data discussed above are read into the program first. As part of this 
process, the maximum rotational quantum numbers are assessed using the algorithm developed 
by Whiting. (See sec. B.2.) These values do not change in the program. 

5.1.4 ■ Continuum Spectrum 

The spectral location of all radiating transitions is determined from the set of energy levels 
in the atoms and molecules appearing in the gas. This information can be used to select a 
minimum set of spectral points for use in the radiation calculation. 

Consider the bound-free continuum. The absorption coefficient for this mechanism is 
a smooth function of frequency with discontinuities known as photoionization edges that 
correspond to the activation of additional energy levels. To resolve this spectrum, spectral points 
must be placed within a small interval on either side' of each jump, as shown in figure 3. However, 
this interval is much smaller than needed for the smooth portions of the curve. The RAD/EQUIL 
code took advantage of this fact, by allowing the user to input a few spectral points to capture 
the jumps. In the LORAN code, this procedure is automated and incorporated in the radiation 
calculation. For the particular set of atomic species considered in a computation, a tailored 
atomic spectrum is generated to resolve the active discontinuities, as the figure illustrates. This 
tailoring simplifies the inputs and eliminates the chances of omitting a level. It also allows the 
atomic continuum calculation (including the smooth free-free continuum) to be performed on 
the minimum resolving array of spectral points. 

As discussed in section 1, the molecular radiation has been modeled here using the smeared 
band approximation that smoothes the rotational structure into an exponential variation. To 
resolve the resulting spectrum, spectral points must be placed around each vibrational band 
head. A few points are sufficient to capture the smoothed rotational variation within each such 
band. The resulting structure is similar to the bound-free spectrum, although it is much more 
complicated because the vibrational bands are far more numerous and often overlap. These 
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Figure 4. Example of spectrum optimization for vibrational band of electronic transition. 

considerations have been used to develop an algorithm to choose an optimal molecular spectrum 
for each input spectral range and set of active species as follows: 

1. For each active species, the vibrational band heads in the spectral range studied are computed 
from equation (88) with J — 0, 

2 . An initial spectrum is generated by resolving the band head and spacing nv points in the 
interval between J = 0 and the maximum rotational level ( J — J rn ax)- An example for a 
major vibrational band of is given by the initial spectrum in figure 4. Note that the 
clustering of points along this curve results from additional weaker bands which underlie this 
major band. 

3. The absorption coefficient decreases rapidly for a given vibrational band (note the log scale 
in fig. 4) so this spectrum can be reduced. This reduction is done by discarding points that 
describe a vibrational band away from its band head (points required to resolve the jump are 
kept), if they are within a given interval crit (a FORTRAN variable) of another point. This 
procedure can reduce the spectral array by nearly an order of magnitude. 

4. A second reduction can be made by discarding any point with a neighbor within a smaller 
interval crit2. 

Application of these steps to the spectrum in figure 4 results in the optimized spectrum shown. 
Application to the entire set of vibrational bands permits a spectrum that resolves the molecular 
transitions with less than 5000 points and that is optimized for any species and spectral range. 
This calculation sifts through many possible spectral points to produce the optimized molecular 
spectrum. It is computationally expensive so it is exercised only once per run. The resulting 
spectrum is stored for all further computations. Appropriate selection of the parameters nv, 
crit , and crit2 controls the amount of detail included in the molecular spectrum and its size. 

After a complete optimized spectrum is obtained, the index of the first and last spectral 
points for each vibrational band is determined. These indexes simplify the logic required in 
computing the molecular band radiation and improve the calculation efficiency. 


17 



About 50 points from the optimized molecular spectrum that resolve strong vibrational 
transitions are identified and interleaved into the atomic continuum spectrum discussed above. 
The resulting spectrum of about 150 points is referred to as the modified atomic spectrum and 
is discussed in section 5 . 3 . 

5.2. Radiation Properties 

5.2.1. Bound- Bound Transitions 

The oscillator strengths f lu , line centers u CL , and Stark HWHM’s 7 ^ are needed to compute 
the bound-bound transitions. The National Bureau of Standards (NBS) has published many line 
centers and oscillator strengths (ref. 27); the Stark HWHM information is available from Griem 
(ref. 3) among others. In addition, the energy levels and degeneracies of the upper and lower 
transition levels are in the NBS compilation. The set of inputs currently used in the LORAN 
code is modified from the atomic line inputs of Park’s NEQAIR code. In addition, for those lines 
for which the HWHM is not known, an approximate correlation by Park is used. (See ref. 28.) 

In the RAD/EQUIL code, Nicolet approximated the so-called high-series lines by an integral. 
These lines have lower state quantum numbers of four or greater, are numerous, and are closely 
spaced below 3.2 eV. Park treats each of these lines individually in the NEQAIR code. In 
the present method, to reduce the number of lines computed and therefore the run time, the 
high-series lines are represented by multiplet-averaged values. Data for these averaged lines are 
obtained from the NBS tables for transitions corresponding to Park’s inputs. The line width 
of each averaged line is then multiplied by an appropriate factor to account for the nearly 
identical closely spaced lines it represents. Integration under this larger wavelength gives the 
total emission and absorption for the multiplet, but transport is calculated for only one line. 
This approximation should have minimal impact on the results, yet it halves the number of lines 
that must be treated. 

The line shape models for the LORAN code are adapted from the development by Nicolet 
(ref. 2 ) as implemented in the RAD/EQUIL code. (See ref. 21 .) Each line is resolved by a 
small number nppl (generally less than or about 15 to 21 ) of spectral points that are distributed 
starting from the line center. (See fig. 1 .) The location of the spectral point farthest from the 
line center is determined from the broadened width of each atomic line and the distance to the 
nearest neighbor. Line-broadening mechanisms treated include Doppler, resonance, and Stark, 
as discussed in section 1 . Lorentz and Doppler line shapes are both computed and the larger of 
b L and is used at each spectral point in a line. If required, additional accuracy is obtained 
from an approximation to the Voigt line shape, combining Lorentz and Doppler effects, which is 
accurate to about 5 percent. 

In Nicolct’s implementation, the line width used to distribute the nppl points was an 
approximation to the maximum width in the gas layer so that all or most of the energy in 
each line would be included. This approximate maximum width was calculated at a location 
input by the user at which the maximum temperature was expected. For an application in which 
a complete flow field is to be computed, the identification of a maximum line width presents 
some difficulties. Further, because the line width may vary by orders of magnitude between 
various points in a flow field, the use of a spectrum geared to the largest width can result in 
significant loss of accuracy. For these reasons, the routine which sets up the line frequency 
spectrum (subroutine freq of RAD/EQUIL) has been modified for application to every point 
in the flow field. This modification results in the specification of a local set of spectral points 
based on the local line width. These disparate frequency spectra are then reconciled to generate 
a single atomic line spectrum. This process results in the selection of a set of spectral points 
capturing the energy in a line at its broadest but that involves detailed spectral information 
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at each grid point. This technique allows satisfactory resolution of the atomic lines without 
requiring excessive storage or computation. 

An additional change made to Nicolet’s method arises because Stark broadening is not always 
the dominant effect on line width. In the original RAD/EQUIL code, the point distribution logic 
in subroutine freq considered only the Stark line width. In the present model, the larger of the 
Doppler or Lorentz (Stark-effect) line width is used. Although this change is small, it should 
result in more complete coverage of the atomic line radiation in the present model. 

The RAD/EQUIL code has also been modified to use nonequilibrium excited state popula- 
tions and electron temperature, where appropriate. The Doppler line width remains a function 
of heavy-particle temperature because it arises from the thermal motion of the atoms themselves. 


5.2.2. Bound- Free Transitions 

The bound-free photoionization edge structure, which is the characteristic shape for this 
mechanism, is produced by the activation of each level in the summation term of equation (13). 
This summation requires determining the lowest accessible bound energy level n* for each 
spectral energy point in the atomic continuum spectrum and summing the radiation properties 
for n* and all levels above it to the ionization limit. The ground-state cross sections are invoked 
as constants to avoid the errors arising from the hydrogenic model, as discussed in section 1. 
Each atomic energy level is treated individually in the bound- free calculation but, because these 
are grouped levels, many more energy levels are in fact approximated and included. Gaunt 
factors are applied to the computation by table lookup of the Biberman and Norman correction 
factor. (See ref. 10.) 

5.2.3. Free- Free Transitions 

The free-free contribution is computed on the same spectrum as the bound-free radiation. The 
resolution thus obtained is more than sufficient and allows the two atomic continuum processes 
to be easily combined. However, because of its form (eqs. (19) and (20)), free- free radiation is 
significant only at low energies. 

5.2.4- Molecular Transitions 

To predict the molecular band radiation, spectroscopic constants must be provided that 
describe the distribution of rotational and vibrational energy levels. These constants are in 
numerous sources, such as Herzberg (ref. 13); Bond, Watson, and Welch (ref. 6); and many 
others. These sources also provide the necessary information on the electronic energy levels and 
energy-level degeneracies in molecules. The intensity of radiation depends on the FVanck-Condon 
factors q {U i and on the electronic transition moments D e Data on these quantities may be 
found in many, often contradictory sources. The complete set of molecular inputs collected by 
Park for his NEQAIR code (ref. 1) has been adopted here for simplicity. 

The molecular band radiation is computed with the optimized molecular spectrum the 
determination of which was discussed in section 5.1. The expressions for the emission and 
absorption coefficients developed in section 4.1 have many terms in common. Advantage is taken 
of these common terms to reduce the computation time. The molecular calculation remains one 
of the more expensive parts of the radiation calculation, so care was taken in programming this 
particular mechanism to ensure that the advantages of vectorization are exploited to the extent 
possible. 
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5.3. Radiative Transport 


The radiative transport calculation employs an approximate treatment of atomic lines based 
on the RAD/EQUIL method. A single mean value of the continuum (atomic continuum plus 
molecular band) absorption and emission coefficient is computed for each line and added to the 
detailed line-only absorption and emission coefficients calculated at the nppl points that describe 
the line. These mean values are obtained by interpolating the information in the optimized 
continuum spectra, thus avoiding an expensive computation at each of the nppl spectral points 
that describe a line. Each atomic line covers only a very small part of the spectrum, so this 
approximation should be fairly accurate. Radiative transport is then computed at all nppl 
spectral points. This detailed result can be integrated over frequency to provide an average flux 
for each line if a concise presentation is desired, or it can be presented in full detail. 

The molecular spectrum can be averaged to obtain values at each point in the much smaller 
modified atomic continuum spectrum, so that the number of spectral points for which transport 
must be calculated is significantly reduced. Because these approximate results are based on 
a detailed molecular spectrum of both emission and absorption, the prediction of molecular 
radiation transport should be quite accurate. 

5.3.1. Numerical Solution 


The integrals in equations (36) -(38) require some care. Because k' v can be negative in some 
regions, the log-linear variation assumed by Nicolet (ref. 2) in the RAD/EQUIL code cannot be 
used. Examination of k' v profiles suggests that a piecewise linear approximation is acceptable. 
For the optical variable, then, the integral is evaluated on the discrete grid as follows: 



For this simple integral, a piecewise linear variation of k ' v is easily shown to lead to the following 
result: 

j n— 1 

T v(Zn) ~ 7^ [ K 'v( z k) + K 'v{ z kJr\)]( z k+\ ~ z k) (42) 

Ar= 1 

For the intensity integrals, the flow-field grid spacing may be such that for strongly absorbing 
spectral regions, A 2 is large with respect to variations in the integrand, particularly because 
of the exponential terms. To integrate accurately, examination of t v profiles suggests that a 
piecewise linear approximation again should be used. By setting ^(s) = a k s + b k in each grid 
interval and by approximating j* based on the average of its values at the two ends of the 
interval, the integral over the exponential term can be carried out analytically. (The obvious 
approach of also setting jf, = c k s + d k causes numerical problems for low values of r„ resulting 
from terms which go to 0/0.) On the discrete grid, the integrals are developed as follows: 


I+(z n ,l) = I up (T w )e- T ^ 


* 1 f z k± 1 / 

+ J2 nlti( z k) + jU z k+l)}e~ Tu(zn] / e a k* +b k ds' (43) 
Jt=l Jz k 

fc max~ l | r z k+\ 

I~(z n ,- 1)= l)]e T " ( * n) / e- a k s '- b kds ' (44) 

U Z J Zu 
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The two integrals remaining in this equation are of a similar form and can be evaluated for 
general and B ^ as follows: 


i: 


k+i A k*+ [i k dx = e n k r +] ( A*dx 

J z k 




n , 


e A k x 


z k + 1 
Z k 


— - J e A k z k+ 1 - 

A k v f 

— — ( e A k z k + \ + B k. - e A k z k + B k\ 

Substituting this result into the integrals in equations (43) and (44), we get 


n — 1 


(45) 


z (Ifc 


k= 1 


A’ max ~ 1 


"max * ^ i v 

-( Z7 „-l)= £ l^k) + 3l^k+l)y vM — (47) 


A' — n 

This expression can be rephrased in terms of the optical variable as 

71 — 1 


I+(z 7l , 1) = I vp {T w )e- Tv{Zn) + J2 ^ [#(**) +#(**+ i)K 


-r t A-n) 


l) - p r id^-) 


A-l 


A'nifix 1 


k=n 


/„ («„,-!) = -]T + j£(**+l)K , ' ( *" ) fj T ' /(; *' +l) - e Tv{zk) 


(48) 


(49) 


The value of n k can be determined by solution of the linear equation t,,(.s) — n k s + b k at the 
two endpoints of each interval as follows: 


T v( z k+\) ~ n k z k+\ + b k 
r c ( z k ) = «k z k + b k 

Subtracting the second equation from the first, to eliminate b k . we have 

T-|/(^-+l) “ T i'( z k) = a k z k + l - (l k z k 
= OA-Cah 1 ~ z k) 

So 


(50) 

(51) 


o-k 


M*k+ 1) - T i'( z k) 


z k+l - z k 

which can be substituted in equations (48) and (49) to get 


(52) 


(53) 
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CUn, 1) - w-" w + E KM +?Az M )} 

e T i'( z k+l) _ e 7/'(**)] 


0 — r i/(zn) 


x e 


^max 1 


( 54 ) 


I-(z n ,~ 1) = ~J2 


1 


^ k -\- 1 % k 


\fu( z k) + jl{ z k+\)] 


k—n 
Tv{zn) 


(55) 


2 [Mn-+ i ) - Mz k )] 

X e — T i/(2fc+i) _ e ~ T i/( 2 Jt) 

This piecewise linear treatment provides excellent results. 

5.3.2. Special Case for Constant r u 

In some spectral and spatial regions, the optical variable may be essentially constant from 
one grid point to the next. When this situation occurs, the numerical integration developed 
above must be revisited. As seen in equations (54) and (55), division by zero occurs in this case. 
The integral in equation (45) can be more simply evaluated for constant Ti/ as 


f Zk +' e A k x+B ki dx w f ZM e A kZk+ B kdx 
Jz k Jz k 

= e A t^B t P + ‘ dx 
J z k 

= e A k z k +B k ( Zk+i _ Zk j 


(56) 


After substituting this result in equations (43) and (44), we get 

w-1 1 

(*«’ ! ) = Up(T w )e~ r ‘'^ + ^ 2 ^ z k) + fvi z k+i)]e Tu(zk) - T ^ Zn \z k+l - z k ) (57) 


k= 1 


1 


^max~A 

/- (z n , -1) = + ti(zk+l)V v{zn) ~ Tv(zk) {zk+i - Zk) 


k~n 


(58) 


5.4. Flow-Field Coupling 

The LAURA algorithm for solution of the Navier Stokes equations is summarized in 
reference 22. The divergence of the radiative flux can be treated as a source term in the two 
energy equations. In general, source terms are handled implicitly (ref. 29) in the LAURA 
algorithm. This treatment improves scheme stability and robustness. Radiation is sensitive to 
flow-field properties, so the added radiation source term should be treated implicitly. However, 
Park and Milos (ref. 30) have found this implicit treatment to be prohibitively expensive in 
computer time except for one-dimensional flows. Furthermore, the necessary terms for the 
Jacobian are nearly impossible to obtain in the current nonequilibrium treatment. Therefore, 
the radiation term is treated as an explicit (ref. 29) source term in the LAURA algorithm. 
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The convergence of the flow field in the LAURA code is characterized by transient waves in 
the gas properties. When the sensitive radiation calculation is performed on a flow field that 
still contains transients, the radiation term can change dramatically. This change can introduce 
additional fluctuations and lead to instability. In the present study, therefore, under '-relaxation 
is applied to the explicit radiation source term. That is, the latest, radiation solution is combined 
with the previous solution to compute an updated radiation term. 

(V • qit) n,+1 = d<i(v • qitf + (i - 4/)(v • qn)" +1 ( 59 ) 

The underrelaxation approach of equation (59) reduces the fluctuations due to coupling and 
increases the stability of the combined algorithm. A relaxation factor d q of 0.5 has been used. 

Because the flow properties change little in a single flow-field iteration and because flow- 
field transients should be allowed to decay, the radiative flux divergence term need not be 
updated with every flow field update. In fact, the total cost of the combined computation can 
be reduced several orders of magnitude by updating the radiation term a minimum number 
of times, as is successfully done in equilibrium shock-fitting methods. Updating it too rarely, 
however, can decrease the convergence rate, produce oscillatory nonconvergence, or destabilize 
the solution. For example, if the flow field changes dramatically between radiation term updates, 
a large change in the radiation will result and can introduce severe waves into the flow field. In 
addition, if flow structures change position significantly between updates of the radiation, the 
radiation term will be temporarily unaligned with t lie associated flow structure. Experience with 
cases characterized by a moderate degree of nonequilibrium suggests that updating the radiation 
calculation once every 500 flow-field updates provides reasonable convergence. For cases with 
more severe nonequilibrium in which coupling has a large effect on the standoff distance, more 
frequent updates are required (every 200 iterations). As a flow field approaches equilibrium, 
on the other hand, the current LAURA algorithm converges more slowly. For near-equilibrium 
cases, the radiation is updated only every 1000 to 5000 flow-field iterations. However, the optimal 
update interval may change if the LORAN code and a different CFD code are coupled. 

Coupled solutions may be difficult to obtain when extreme gradients occur in the shock wave, 
as is the case for near-equilibrium flows. The computational mesh used in a flow-field calculation 
seldom provides adequate resolution in the shock region to produce smooth temperature profiles 
in such flows. This lack of resolution at the shock does not significantly degrade the flow-field 
solution. The radiation calculation is sensitive to temperature, however, so it can cause large' 
changes in the profile of V • RlL which may lead to inaccurate results or even divergent solutions. 
Increasing the number of grid points can alleviate this problem but is not always practical 
because of storage limitations. Grid adaption, on the other hand, can improve the temperature 
profiles without increased storage. To provide this ability, a grid adaption algorithm has been 
added to the LAURA code. (See ref. 31.) This algorithm can be used to cluster grid points 
around a specified density or temperature in the shock. 

The low-temperature (precursor) region ahead of the shock in a shock-capturing scheme also 
requires special treatment. The LORAN radiation model as currently implemented does not 
include accurate radiation properties at low temperatures. The precursor is generally assumed 
to have negligible effect on the converged flow field as it can absorb only a small fraction of 
the radiative energy passing through it. Recent results (ref. 32) confirm the validity of this 
approximation. Therefore, V ' RR '° zero the value appropriate for a nonparticipating 

medium when the vibrational or translational temperature is below a certain minimum. This 
adjustment is done within the LAURA flow-field calculation every iteration. Even if accurate 
precursor (low-temperature) radiation properties are available, some logic is required at each 
iteration so that a decrease in shock standoff distance between radiation updates will not leave 
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a large radiation term in what is now the free stream. A temperature limit of 500 K provides 
good results for the cases so far examined. 

As a practical matter, a converged flow-field solution should be obtained first with the LAURA 
code, then the radiation should be “turned on” to converge the coupled flow field. If the radiation 
is only a small perturbation to the flow field, this second convergence may require only a few 
additional LAURA iterations. If strong radiation effects are evident, then the perturbed flow 
field will require many additional iterations to converge. The optimal sequence for converging a 
coupled solution is not known; it must be learned by experience with the code. 

6. Computational Optimization 

Obtaining radiation predictions for complete flow fields or coupling radiation and flow-field 
solutions requires that the radiation calculations be carried out many times. To make this 
iteration practical and to minimize the cost, the radiation program should be optimized and 
the minimum set of calculations which will provide an accurate radiation prediction should be 
determined. The accuracy required will vary depending on the flow regime and body shape 
and on the intended use of the results. Therefore, an optimization study may be required for 
each project undertaken with the LORAN code. (An example is given in ref. 33.) Methods of 
achieving this optimization are discussed in sections 6.1. 6.3. 

6.1. Radiation Calculation 

Parametric studies can be conducted to determine the spectral resolution required for 
acceptable accuracy in the predicted radiation. The effects of the following variables should 
be studied: na, the number of points allowed in the atomic continuum spectrum; nppl\ and nu, 
crH, and crit2 , the several variables that control the generation of the spectrum for molecular 
band radiation. In addition, various sets of atomic lines can be inserted, from the fully detailed 
to a set of multiplet-averaged lines. The parameters and inputs should be chosen as appropriate 
for each study. 

6.2. Excitation Calculation 


The computing time required to obtain the nonequilibrium excited state populations is a 
significant fraction of the total for the LORAN code. The QSS algorithm used to perform this 
calculation was obtained from the NEQAIR program and is used as a “black box.” Improvements 
to this algorithm could significantly improve the computing efficiency of the radiation prediction. 
However, no changes have been made to date to take advantage of that savings. 

6.3. Radiation Subgrid 


Calculating radiative properties at every point on a flow-field grid is often wasteful. The 
absorption and emission coefficients change significantly only in regions with large gradients in 
temperature or species concentrations. Bolz (ref. 34) developed an algorithm to automatically 
select a subset of grid points for the radiation calculation. He defines a weighting function Z 
which can be adapted to the nonequilibrium situation as 
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Tj , cm 

Figure 5. Example of radiation subgrid selection. 

The prime in equation (60) denotes the partial derivative with respect to the normal to the 
wall. This quantity is evaluated at each normal grid point k in the flow field. If the vibrational 
temperatures of individual species are computed in the future, additional terms can be added 
to Z k . The average derivative of the species number densities is 


AT' = - 

Tic 


n s 




3 - 1 


(61) 


for s chemical species. Radiation is then computed for a subset k\ { of the flow-field grid points 
chosen at equal intervals of Z k . An example of such a grid is given in figure 5. Although 
weighting factors W, are included for each term in Z k , the gradient in one variable should not 
necessarily be emphasized more than that of another. Therefore, these factors are all set equal 
to 1. In some applications, different values of the weighting factors may benefit the user. For 
optimal accuracy in a study, the size of the radiation subgrid k\ { should also be varied. (See 
ref. 33 for an example.) 

To further reduce the required computational time, every other normal line of the flow-field 
grid and every other plane in a three-dimensional solution might be skipped in the radiation 
calculation. The required radiation properties can be interpolated for the skipped normal lines. 
Although the subgrid algorithm could be applied in these directions also, the expected gradients 
are gcnerallv gentler and not as well suited to the method. 


7. Using the Code 

The code consists of 58 FORTRAN subroutines and 6 “include files” and has been run on 
Sun, Silicon Graphics, and Cray computers, all with UNIX operating systems. A “makefile” 
is available for compilation and loading. Double precision accuracy is required on the Sun 
and Silicon Graphics computers; on the Cray computer, single precision is adequate and much 
faster. Conversion to single precision requires removing implicit double precision statements 
everywhere, changing double precision declarations to real, and switching all intrinsic functions 
to a single precision form; otherwise, the user must employ an appropriate compiler option if 
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available. For a brief overview of the code, appendix C presents the code structure in the form 
of a subroutine calling tree. 

The code is provided to potential users with sample input and output files to allow verification. 
Because these files are quite extensive, they are not included in this document. Examples of 
results obtained with the LORAN code can also be found in references 18, 26, 31, 33, and 35-38. 

7.1. Include Files 


The six include files set several important parameters for the code and have comments with 
descriptions of each variable. The following are the general functions of the include files: 


btask . zzzz 
geoa. zzzz 
parml .zzzz 

parm2 . zzzz 
parm4 . zzzz 
parr . zzzz 


storage for heat flux and flux divergence results 
storage of flow-field grid geometry 

parameters that set size of flow-field grid; must be changed for different 
grid sizes 

parameters that set number of species; no change if air 
grid singularity flag 

parameters that set radiation subgrid size and control optimized radia- 
tion spectrum; change for different grid sizes and spectral resolution 


7.2. Input Files 


Two input files are required: 

neq.in contains excitation and radiation data adapted from NEQAIR input files 

radneq. in contains control and line data adapted from RAD/EQUIL input 

These files are opened in the code so need only be present in the directory where the code is being 
run. An input guide is included in appendix D. Except for a few variables near the beginning 
of file radneq. in, most of these variables will not be changed by the user. The description of 
the variables is provided to aid in understanding the code and for adding to or updating the 
excitation and radiation properties. 

Although the input has provision for 25 species, with 22 ordered as in the RAD/EQUIL 
code and with 3 spares, the LORAN code does not yet contain radiative properties for all these 
species. Some properties are available for N, O, C, N + , 0+, C+, e“, N 2 , 0 2 , NO, C 2 , CO, 
CN, and N^. For N, O, and C, bound-free and line radiation is computed. For N + , 0 + , and 
C + , free- free radiation is considered. The following molecular radiative transitions are included: 
); N 2 (l + ); N 2 (2 + ); N 2 Birge-Hopfield; NO /?, NO 7, and 0 2 Schumann-Runge; CN violet; 
CN red; CO(4+); and C 2 Swings bands. As more properties become available, they will be added 
to the LORAN input set. 

Flow-field condition input is also required through some interface. In the stand-alone LORAN 
code, these data reside in the file restart.dat, which is also opened in the code. 

7.3. Output Files 


The following output files may be produced, depending on the control flags set: 

standard output name may be set by user; contains FORTRAN errors and warn- 
ings; echoes parameter values used; summarizes spectrum and flux 
calculation 
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flow. out 


fluxdiv.out 


qr . out 


rademis . out 


single-zone I-ordered Tecplot 2 computer file that contains normal 
coordinate and flow-field properties along stagnation grid line with 
4 temperatures and 11 species number densities (for air) 

single-zone IJ-ordered Tecplot file that contains grid coordinates and 
divergence of radiative flux for each normal grid line 

multizone I-ordered Tecplot file that contains normal grid coordinate 
and positive and negative intensity /flux ( // and I u or '/,! and q u ) for 
each normal grid line (one zone per grid line); uses block format 

single-zone IJ-ordered Tecplot file that contains grid coordinates and 
radiative emission at each point in radiation subgrid; includes atomic 
continuum, molecular, atomic line, visible, and infrared-only atomic line 
emission 


radneq. out 


walrad . out 


summary of data input and radiation calculation; gives total emission 
results and reports irregularities in calculation 

single-zone I-ordered Tecplot file that contains path length coordinate, 
upper limit of spectral range, and wall radiative intensity/flux at each 
surface grid point 


wsp.out multizone I-ordered Tecplot file that contains spectral output. 

frequency points, spectral intensity/flux, and running integral of 
intensity/flux for each normal grid line (one zone per grid line) 

Tecplot files are the input files for the Tecplot plotting program. (See ref. 39.) A description of 
their format is given in appendix E. 


8. Conclusions 

A new method for computing radiative heating in thermochemical nonequilibrium flow fields 
has been developed. This method takes advantage of a priori knowledge of the radiation spectrum 
to generate an optimized spectral array for radiation calculations. Computational efficiency is 
further enhanced by the assumption of smeared molecular bands. The theoretical basis and 
practical application of the method have been summarized here. 


NASA Langley Research Center 
Hampton, VA 23681-0001 
May 23, 1994 


2 Tecplot: trademark of Amtec Engineering, Inc., Bellevue, WA, 98009-3633. 
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Appendix A 
Gaussian Units 

For the reader’s convenience, the Gaussian system of units is redefined here. 

c speed of light in vacuum, 2.9979 x 10 10 cm/sec 

e electron charge, 4.80286 x 10“ 10 statcoulomb 

h Planck’s constant, 6.6262 x 10 -27 erg-sec 

k Boltzmann’s constant, 1.3807 x 10“ 16 erg/K 

m e electron mass, 9.1095 x 10“ 28 g 

The following relations hold among several of these quantities: 

The Bohr radius «o is given by 

h? 

0 47r 2 m e e 2 
= 0.52918 x 10 _8 cm 

The fine structure constant a, a dimensionless quantity, is 

27re 2 

Q ~~ ~h^T 
= 7.292 x 10“ 3 
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Appendix B 

Molecular Radiation Model 


Smeared Band Derivation 


In the discussion that follows, Zeldovich-Raizer refers to reference 9; Herzberg refers to 
reference 13. The energy differentials for the molecular transitions are established first. The 
expressions given below for the energy levels are to some extent based on curve fits. In most 
cases, improvements are possible. 

For electronic transitions, the energy levels are given by E e , the electron energy of a particular 
energy level measured from the ground state. For vibrational transitions, the energy levels are 
obtained from the fourth-order expression given by Herzberg, 


E 1} = he 


uJ e ( V + ~ UeXe + +Weye ( v+ 2) +UJeZe (^ v+ 2 ) 


(Bl) 


which includes a correction for the anharmonicity of the vibrational energy potential. The 
vibrational levels range from the zero-point energy at v = 0 to some maximum v determined 
by the dissociation of the molecule. The coefficients in the above expression are spectroscopic 
constants and are available from several sources. The values used in the LORAN code are 
those collected by Park for the NEQAIR code. The energy levels for the rotational transitions, 
including a correction for nonrigid oscillators, are 


E r — he 


B V J(J + 1) - D V J 2 (J + l) 2 


(B2) 


where the spectroscopic constant B v has been corrected for rotational effects as 



and the second rotational spectroscopic constant D v is 



These levels range from zero energy at J = 0 to some maximum J, again determined by 
dissociation of the molecule for each vibrational level. The spectroscopic constants for rotational 
energy states are also those collected by Park. 

The complete molecular energy level is 

E tot — E e -f E v T E r (B5) 


where the individual energy contributions are given by equations (Bl) and (B2) above. Desig- 
nating the upper level by the superscript u and the lower level by l, the wave number (inverse 
wavelength) of a particular transition is 

1 _ £?ot ~ E L (B6) 

A he 
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For convenience in manipulating these expressions, one special wave number is defined, the wave 
number with zero rotational quantum number and arbitrary vibrational quantum numbers. 
Thus, 

1 - E l ( , El 1 - El 

t = ~^r + ^~r (B7) 

\u v i he he 

The formidable expression for the wave number of an individual transition is simplified by 
the selection rules for rotation. (See refs. 9 and 13.) A detailed discussion of these selection rules 
and of the designation of molecular states is outside the scope of the present work. Nevertheless, 
these rules allow the development of three specific wave numbers for the transitions A J = — 1, 
0, and 1, which are denoted P, Q, and R, respectively. Therefore, 

jh = + ( J 1 ) A (d[ - D^+{J l f(2D[ + 2D“) 

+ (J 1 ) 2 (d[ - B[ - £>“ + 5«) + J l (-fl“ - B l v ) (B8) 

~r~ = T-h- + {J l ) A (D l e - D»)+(J l f(2Dl - 2 D?) 

X Q \,u v l V > V ) 

+ (J 1 ) 2 (dI - B l v - Dg + B “ ) + J 1 (b;; - B l v ) (B9) 

T = + {J 1 ) a (d[ - D^+(J l f(2D l e - 6 Z?“) 

+ (J 1 ) 2 (d[ - b{, - 13 D“ + B^ + J 1 ( 3 - 12/J>" - + 2Z?;; - 4D f M (BIO) 

After the molecular band spectra are thus determined, the probability of a transition between 
levels must be evaluated to obtain an expression for the absorption coefficient. This probability 
is obtained from a quantum mechanical analysis. The resulting expression for the transition 
probability is given in Zeldovich-Raizer: 


A ul 


= A 


Bv u J u 

Av‘j‘ 


64tt 4 / \3 2 

= 3gnhct \Bv u J u ,Av l J l ) D c\BA%n v iPjuji 


(Bll) 


where the upper state degeneracy g $ has been added to the denominator to account for 
degenerate states. (See ref. 13, p. 21.) Also, q vUyl is the Franck-Condon factor for the probability 
of a vibrational transition, p JU ji is the probability of a rotational transition, and D^ B . is the 
square of the matrix element (sre 2, used by Park in the NEQAIR code, is the nondimensional 
value, which is divided by (aoe) 2 , the squared product of the Bohr radius and the charge of 
an electron) arising from the wave function in the quantum mechanical analysis. The Einstein 
coefficient for spontaneous emission is A and super- and subscripts A and B denote the lower and 
upper electronic states of the molecule, respectively. Using the principle of detailed balancing 
leads to the well-known relationship between the Einstein coefficients for emission and absorption 
(i.e., refs. 27 or 9), 


Bln - 


C 1 ju A 

hu l.l 91 


(B12) 


If this relationship is expressed for energy states rather than energy levels, the degeneracies g u 
and gi are identically one, by definition. As a practical matter, the use of energy levels, which 
may be degenerate, is preferred because it reduces the bookkeeping required. 


30 



The absorption coefficient k/ u is related to the Einstein coefficient for absorption B iu (ref. 40) 

by 

R’lu niB iu hisfoi 

Substituting equations (Bll) and (B12) in this last expression yields the desired expression for 
the spectral absorption coefficient, 

^AvW.BvVJv = 9 g^J jt ^^ N Av‘j^Avlj‘,Bv-J^ D lBA%n v lPjujlF(^ (B14) 

where F{v} is a normalized line shape. This expression must be evaluated for each .4 ■ B , v , 
v l , J u , and J l to obtain the complete absorption spectrum of a molecular species. The smeared 
band model averages out the rotational line structure by summing equation (B14) over J u and 
introducing an average frequency V ByU Av i for the rotational lines. The result is 
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U Av l J l ,Bv u 


8 tt 3 J_ 
3 he <74 


N Av^ jl BA9 v u v l ^ ^ Bv 11 ,Ai/ 


(B15) 


where use has been made of the relation Pj u jl — 1> which expresses the fact that all 
transitions occur between two levels of the molecule. The degeneracy has been expanded into 
its components, 

9AvJ = 9A9v9j ( B16 ) 

where 

9v = 1 (B17) 

because no degenerate vibrational levels exist and 

<7j = 2J + 1 (B18) 


For most rotational lines, J > 1 and the ratio of rotational degeneracies in equation (B14) is 
approximately equal to 1 . This approximation has been used even for transitions with low values 
of J. The error introduced in a single transition with a low rotational quantum number can be 
as much as 300 percent; however, only a few transitions have low J, so that the effect on the 
overall absorption coefficient is small. The error incurred by introducing the average frequency 

„ , / for all rotational lines in the band v l - v u is minimal, to 7 percent. 

Dv u ,Av l 

To obtain the complete absorption coefficient, equation (B15) is now summed over the final 
vibrational level v u . This summation yields 


87T 


V Av l J l ,B 


She <74 N Av l J lD l\BA F ^)'^2 / %^v ll 'Bv u ,Av l 


(B 19) 


Here, the frequency remains inside the summation and is not averaged, as the variation can be 
large. 

The mean value theorem is then used to define an average absorption coefficient for a line by 
integrating over one line width for 
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-'line 
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~ N Av l J lD< e\BAY^ q v u v lV Bv u ,Av l = Ku Av l J l .B ^ 


(B20) 


where Av is the line width. 
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Equation (B20) is then solved for k u and summed over the initial rotational and vibrational 
levels to obtain 

O 3 1 OO OO 0° *r 

< B2 » 

V* jl 

where the sum term is zero for any transition which does not absorb radiation in the frequency 
interval considered. 

Equation (B 21 ) is the complete formula for the smeared band absorption coefficient for a 
molecular species. In this form it contains two approximations: an average frequency for the 
rotational lines in a vibrational band and the averaging over a single line. To use this expression 
in a radiation model, the various parts must now be amplified in terms of known quantities. 

In nonequilibrium conditions, the number density in a particular energy level is given by the 
following expression, where the existence of rotational and vibrational temperatures and the 
separability of the rotational and vibrational energy modes have been assumed: 


N 


Av‘j‘ = N * 


exp (-E^/kTy) gj 

Qv(T v ) 


ex P (-Eji/kTr) 

QATr) 


(B 22 ) 


where Q v and Q r are the vibrational and rotational partition functions, respectively. The 
population of the electronic energy level /V 4 must be obtained from some nonequilibrium 
excitation calculation. In the present work, Park’s QSS method (ref. 14) has been used. The 
partition functions are obtained from their definition 


Qi(Ti) = ^]exp (-Ei/kTi) 

i = 0 


(B23) 


In practice, the upper limit of this summation is imposed either by predissociation in the 
rotational case or by numerical limits arising from inaccuracies in the higher order coefficients 
of equation (Bl) for the vibrational case. (The rotational case is discussed in a separate section 
below.) For the vibrational case, the sum may be terminated when the contribution of the next 
term is less than some threshold percentage of the partition function. The size of each term 
must also be monitored to detect the onset of significant inaccuracies in equation (Bl); these 
inaccuracies occur at a high enough vibrational quantum number that the sum for the partition 
function can be truncated at that point. 

The next term to be examined in equation (B21) is A v, which introduces another approxi- 
mation of the smeared band model. The line width Az^ is assumed to be approximately equal 
to an average line spacing. To evaluate the line spacing, consider the three expressions for the 
line centers, equations (B 8 ) (BIO). If the terms D c are neglected compared with the terms B v 
(they are much smaller), the equations reduce to 
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For a line with J l » 1, which is the case for most lines, these three expressions collapse into the 
single form 


1 v 
A c 

« - 1 — + (J 1 ) 2 ( BJ} - Bl) (B27) 

V U V< V ' 

Solving for v and taking the derivative with respect to J (remembering that A J = 1 for adjacent 
lines) gives 

Av = 2cJ l \B v u - B v i\ (B28) 

where the absolute value is taken to ensure that the average spacing is positive. Equation (B27) 
also yields an approximate expression for J 1 in terms of the wavelength or wave number 
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Because (J^)^ > 0, the right side of this equation must also be positive. This situation implies 
for B v u > B v i that X yUyl > A and for B v u < B yi that A yUy i < A. These requirements are the 
mathematical statements of the fact that each vibrational band absorbs in only a portion of the 
spectrum. In an evaluation of the absorption coefficient for a particular band, the wavelength or 
frequency considered must be checked for the above conditions to avoid superfluous contributions. 

The rotational energy levels in equation (B2) contain the assumption that J* 1. This 
assumption allows the replacements J* - K 1 ~ J 1 and 2J^ + 1 ~ 2 J 1 . Equation (B29) can then 
be substituted in equation (B2) with the result 
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This expression can then be substituted in the rotational exponent in equation (B22). It 
constitutes the final approximation of the smeared band model. 

The last term to evaluate in equation (B21) is the averaged frequency V Bv u Av l ■ This term 
must be evaluated separately for so-called S and non-E transitions according to the rotational 
selection rules mentioned above. For a given vibrational transition, equations (B8) (BIO) 
describe the line centers of each rotational transition (eq. (B9) is not allowed for E transitions). 
Because the 0-0 rotational transition is always forbidden, Vjj v u Av i can be evaluated as follows: 
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Bv u ,Av 


1 ( A ) Bv n ,Av^ 



(B32) 


where equation (B31) is for non-E transitions and equation (B32) is for E transitions. 

Substituting equations (B8) (BIO) in equations (B31) and (B32) and collecting terms result 
in the following expressions for the averaged frequency: 
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where again the first form is for non-E transitions and the second is for E transitions. These 
equations could be implemented as is; however, the result is simpler and still accurate if the 
terms D e are neglected compared with the terms B v . The result is 
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V Bv u ,Av l 
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(B36) 


These equations for non-E and E transitions can then be substituted in equation (B21) and the 
averaged spectral absorption coefficient can be calculated. 
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The final result for the smeared band molecular absorption coefficient is 
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where u is given by equation (B35) for non-S transitions and by equation (B36) for S transitions, 
and the vibrational energy E j is evaluated from equation (Bl) with the appropriate values for 
the lower energy level. To be consistent, the approximation 2 ,7^ T 1 ^ 2.7 is used to reduce the 
fraction inside the summation over .7^. 

Alternatively, the mean value theorem can be applied in equation (B14) to obtain an average 
absorption coefficient for a single rotational line. The same summations and averaging as above 
can then be performed. The only difference from the result in equation (B37) is that the averaged 
frequency v j iv u ,\ r l i s replaced by i ' Av l jl jj v njh which is essentially the speed of light times the 

wave number A appearing in the exponential term above. This substitution eliminates the 
need for the calculation of V RvU Av i ■ The approximations required to obtain this result are that 

v Av l J l ,fiv u J l ~ u Av l J 1 .Bv u J 1 J \ 

~ l 'Av l J l J3v"J>-\ 

which must, be true for the closely spaced rotational levels, and that 
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f or ./" _ jl % jn — jl _ i an d ju - jl + j As discussed above, this approximation will be 
true except for a few low values of ,7 and is acceptable. In general, it will introduce less errors 
than would the calculation of v Bv u This alternate smeared band development is used in the 
LORAN code for both emission and absorption, although the coding for both has been retained. 
The smeared band absorption coefficient is then 
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A smeared band expression can be developed for the molecular emission coefficient by a 
process parallel to that for absorption. The total emission in a rotational line J* d is obtained 
from the Einstein coefficient for spontaneous emission: 


J ul = N uA u ihv ul 


(B41) 


The emission intensity at a particular frequency is obtained by dividing by the solid angle and 
introducing F{v). Thus, 

N u A ul his ul 

-F{v) (B42) 


Ju ul ~ 


4n 


Substituting equation (Bll) for the Einstein coefficient yields 
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To obtain the smeared band result, this expression must first be summed over J l and the average 
frequency i ' Bv u Av i introduced. Using the summation rule for p JUJ i gives 
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In parallel to the development for absorption, the sum over v l is now performed. Thus, 

167T 3 00 

i v A,BvUJ* = N Bv u J u 3 5fiC 2 D 1\BA F (v) ^2 Q v u v l [ l/ Bv u ,A v 1 ) 4 ( B45 ) 

v l 

where again the frequency remains inside this summation for accuracy. The average emission 
coefficient for a single rotational line is now obtained using the mean value theorem. 
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This equation can now be solved for the rotationally averaged emission coefficient. The 
complete smeared band emission coefficient is then obtained by summations over the final 
rotational and vibrational levels J u and v u . The result is 


BA E E w E ^ (B47) 

V U v l JU 

As in the case for absorption, the infinite sums are truncated by dissociation. 

The average frequency v Bv u Av l in equation (B47) is the same as that for absorption and is 
given by equations (B35) and (B36) for non-E and E transitions, respectively. The average line 
spacing is also the same and is given by the expression in equation (B28) with the rotational 
quantum number expressed in terms of the wave number according to equation (B29). The final 
term in the emission coefficient, the nonequilibrium upper level population, is found from the 
expression 

N Bv uju = N B 9ju exp(-Eju/kT r ) 

Qv(T v ) Q r {T r ) (B48) 
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where the vibrational and rotational energies are given as a function of wave number by 
equations (Bl) and (B30), evaluated for the upper energy level, and the electronic level 
population N B is obtained from the QSS method. The evaluation of the partition functions 
was discussed in the section on absorption. 

Substituting all these expressions in the emission coefficient of equation (B47) yields the. 
smeared band result 
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Again, the approximation 2.7“ + 1 » 2.7“ is used inside the summation over J u for consistency 
with the level of accuracy of the result. 

As mentioned above, an alternate development is possible in which the mean value theorem 
is applied to equation (B43). This alternate model in which v Bv u Av i is replaced by l/ Av lji Bv ujl 
is the one implemented in the LORAN code, although coding for both has been retained. The 
smeared band emission is then calculated from 
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Maximum Rotational Quantum Number 

The maximum rotational quantum number J ma x can be determined to varying levels of 
approximation. The first-order calculation would predict the maximum when the energy of a 
level exceeds the dissociation energy of the species. As discussed in Herzberg (ref. 13), however, 
predissociation may occur due to the shape of the energy potential curve. Whiting, Arnold, and 
Lyle (ref. 41) developed a computer code based on this fact to predict the maximum rotational 
quantum number. The method was later refined by E. E. Whiting, Eloret Institute, Sunnyvale, 
California, in a private communication. This most recent method has been adopted for use here. 
For every vibrational band of each molecular contributor, it computes the maximum allowable 
rotational quantum number and is based on the shape of the Morse centrifugal potential that 
describes the rotating and vibrating molecule. 
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Further enhancements to the theory for the maximum rotational quantum number are 
possible, but this level of accuracy is sufficient for the current purpose. 
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Appendix C 
Code Structure 


Flowtrace Calling Tree Report 


Ordinal 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 


Level 

1 

2 

3 

4 

4 

5 

6 

5 

6 
7 
6 

5 

6 

7 

8 
6 
7 
7 

7 

8 

7 

8 
9 

10 

8 

8 

9 

10 

9 

9 

9 

8 

8 

8 

6 

7 

7 

5 

6 
6 

5 

6 
6 
6 
5 


Routine name 
+SYSTEM 
RMAIN 
RADIAT 
EVENLN 
RDN0NEQ 
ASPECT 
SORT 
C0NTN2 
MUC0NT 
TAINT 
MUM0L 
EXCT 

MDISSEQ 

EQC 

PARTM 

MOLEXT 

EQC (REPEAT to Tree Ordinal 14) 

PARTM 

PARTQ 

PARTM 

QSSM 

ELECD 

CROSAB 

TAINT 

EQC (REPEAT to Tree Ordinal 14) 
EXCTML 
AI3 
SIMP 

CROSAB (REPEAT to Tree Ordinal 23) 
GAMC 
GAMI 
EXCTMM 
PMINV 
RATED 

QSS 

EONE 

PMINV 

LINT2 

FREQ 

MULE 

MSPECT 

CULL2 

CULLER 

MSORT 

NODEN 
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46 6 

47 5 

48 6 

49 7 

50 8 

51 7 

52 8 

53 5 

54 4 

55 5 

56 6 

57 4 

58 5 


PARTM 

RADIN 

MOLIN 

JMAX 

FKMAX 

VMAX 

FKMAX 

SORT 

SLABTR 

TSFULL 

TRNEQ 

SUBGRD 

TPD 
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Appendix D 
Input Guide 

Flow-Field, Control, and Line Inputs 


Flow-Field Data : TV onuni f orm Condition Inputs 

Below are the inputs for the flow field. This group contains the nonuniform condition inputs. 

The following information is passed through the call statement to the radiation subroutine 
rdnoneq. 

etaln(k) array containing the normal coordinate of the grid points 

on the current normal line of the flow-field grid, cm 


tln(k) , trln(k) ,tvln(k) ,teln(k) arrays containing respectively the heavy-particle trans- 
lational, rotational, vibrational, and electron-electronic 
temperatures along the normal grid line, K 

(xnln(k, j) , j=l ,25) array containing the number densities in cm' 3 of the 

radiating species along the normal grid line in the 
following order: N, O, C, H, Ar, N + , 0 + , C + , N , O , 
C“, H“, e“, N 2 , () 2 , NO, C 2 , H 2 , CO, CN, C 3 , and Nj, 
with three spares; the program can compute the number 
densities of the negative ions N”, 0~, C , and H using 
a Saha equilibrium relation for these trace species; it also 
can compute an equilibrium number density of if no 
value is included for it 


kir 

ispet 


k dimension of radiation subgrid along normal line 

flag for optimized spectrum computation (initially set to 
(); code sets to 1 after spectrum is computed) 


Control Data 

The following data are contained in file radneq.in. They are read into the program in 
subroutine radin.f. Generally, only the values on the first five cards in this input file must 
be changed unless a user has newer atomic line data. (Note that these inputs are called cards; 
in fact, they are simply lines in the input file. The terminology has been retained from the 
RAD/EQUIL code.) 


Card 1, Format (i‘10.4) 

Field 1, (Columns 1-10), tw 

Wall temperature (K). Used in the radiative transport calculation. To 
obtain the incident radiation, this value should be set to zero so that 
the wall reradiation does not affect the result. 
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Card 2, Format (2fl0.2) 

Field 1, (Columns 6-15), wmin 

Minimum frequency for the radiation computation (eV). 
Field 2, (Columns 16-25), wmax 

Maximum frequency for the radiation computation (eV). 

< J 


' v 

Card 3, Format (20il, 15a4) 

Field 1, (Columns 1-20), kr 
see description below 
Field 2, (Columns 21-80), CASE 

Title of the case (alphanumeric) used for identification of printed out- 
put. 

> 


Description of kr flags 
Column 1 not used 


Column 2 
0 
1 


determines whether continuum radiation is included 
include 

do not include 


Column 3 
0 

1 

Column 4 
0 
1 


determines if molecules are to be included in the calculation 
molecules can be included (see also the ncrc variable which overrides 
this control) 

molecules are not included (regardless of what ncrc is) 

determines whether the line contributions are calculated 
exclude lines 
include lines 


Column 5 

0 

1 


determines whether equilibrium or nonequilibrium excited state popula- 
tions are to be used 
nonequilibrium 
Boltzmann equilibrium 


Column 6 
0 

1 


determines accuracy of rotational partition function used 

approximate 

exact 


Column 7 
0 
1 


determines whether intensity or flux is to be calculated 

intensity 

flux 
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Column 8 
0 
1 

Column 9 


0 

1 

Column 10 
0 
1 

Column 11 
0 
1 

Column 12 
0 
1 

Column 13 

0 

1 


determines whether frequency or wavelength is used in output 
uses frequency, his, eV 
uses wavelength, A, A 

determines if the program is to check the validity of the continuum fre- 
quency grid (i.e., do the input spectral limits include at least 95 percent 
of the radiative energy?) 

the grid is checked, use for energy transfer calculations 

the grid is not checked, use when only a selected part of the spectrum is 

of interest 

determines whether incident or transmitted radiation is to be calculated 

incident 

transmitted 

determines whether tangent slab transport is used 
nontangent slab 
tangent slab 

determines atomic line shape used 

Doppler /Lorentz 

Voigt 

determines whether surface absorption and transmission properties are 
to be read in 
not read in 
read in 


Columns 14 to 16 

Column 17 
0 
1 
2 


not used 

determines the amount of output 
normal output 
moderate output 
extensive output 


Columns 18 to 20 not used 


Card 4, Format (24i3) 

Field 1, (Columns 1-3, 4-6, etc.), ncrc 

Determine which radiation phenomena are to be included in the calcu- 
lation and which are not. Flags for the species are given in the order 
shown in Table D.I, with meaning shown at right. Note that for the 
molecular species, the code is currently set up to include all bands or 
none. There is no provision to turn single molecular bands on and off 
from the input file. See the section on Input Files for a list of the 
molecular bands currently in LORAN. 
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Card 5, Format (6el2.8) 

Field 1, (Columns 1-12), tob 

Estimated maximum flowfield temperature (K), used to determine the 
necessary spectral region (see kr(9 ')) LORAN will also use the maximum 
temperature input in the Flowfield Data described previously. 

> 

Table DI. Radiation Phenomena Flags (See Card 4) 


Index 

Species 

Flag for species | 

Atomic and ionic species 

1 

N 


2 

O 


3 

C 


4 

H 

0 — * skip entirely 

5 

Ar 

1 — * bound free 

6 

N+ 

2 — > bound free 4- line 

7 

0+ 

3 — ► bound free T free free 

8 

C+ 

4 — ► bound free T free free + line (every atomic type) 

9 

N" 

5 — ► free free + line 

10 

cr 

6 — ► free free 

11 

cr 

7 — > line 

12 

H- 


Electron 

13 

e j 


Molecular species 

14 

n 2 


15 

O 2 


16 

NO 


17 

C 2 


18 

H 2 

► 0 — > not included 

19 

CO 

1 — * include through subroutine mumol 

20 

CN 

(Smeared rotational line model) 

21 

c 3 


22 

N+ 

2 > 



Line Data 

Basic radiation data. This group is the atomic radiation database and should not be changed 
unless new species are added, new data become available, or a different treatment of multiplet 
lines is to be used. Special treatment is given to the four Lyman and Balmer a and (3 lines. 
Half-width information (a gamp value; see “Data on Line Transitions”) is read in for each of the 
remaining hydrogen lines. (Note that the half-width information for the Lyman 7 and Paschen 
a lines is used to set up the frequency coordinate system only. For the lines which are still higher 
than these, the input HWHM is used throughout the transport calculation.) The hydrogen line 
calculation is from Nicolet’s work and has not been checked for the LORAN code. Therefore, 
the user should expect problems. 
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Card 1, Format (40i2) 

Field 1, (Columns 1-2), nhv 

Number of line groups to be included in the calculation (maximum of 
25). A discussion of the judgments involved in selecting the appropriate 
number of line groups is given in the discussion of the fhvm and fhvp 
values. 

Field 2, (Columns 5-6), nxi 

This is the number of special hydrogen lines for which no HWHM data 
is read in. When the hydrogen lines are to be included, nxi = 4. Oth- 
erwise, it equals zero. 


Line group specification. 


Cards 1, 2...nl, Format (6el2.4) 

Field 1, (Columns 1-12, 13-24, etc.), fhvm 

Low frequency boundaries on the line groups. An fhvm value along 
with an fhvp value defines a frequency increment. All the lines with 
their centers within that increment define a line group. The continuum 
properties are evaluated at only one frequency point for each line in 
the line group. Contributions to the transport from sources outside the 
boundaries of the line group area are not taken into account. Therefore, 
the boundaries of a line group should never be too close to the center 
of a line within the group. 


Cards nlTl, nl+2...n2, Format (6el2.4) 

Field 1, (Columns 1-12, 13-24, etc.), fhvp 

Corresponding high frequency boundaries on the line groups. 
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Cards n2-fl, n2+2...n3, Format (6el2.4) 

Field 1, (Columns 1-12, 13-24, etc.), fhv 

Frequencies used to determine whether a line group falls within wmin 
and wmax. Usually taken to be roughly midway between fhvm and 
fhvp. Thus fhv is referred to as the “average” frequency of the line 
group. 



Hydrogen index ( skip this group if nxi = 0). 


> 

Card 1, Format (40i2) 

Field 1, (Columns 1-2, 3-4, etc.), ia 

Indices on the line groups which contain the special hydrogen lines (one 
per group, maximum). The line groups are always numbered sequen- 
tially starting with those at the lowest frequency. (The calculation of 
hydrogen lines is from RAD/EQUIL and has not been checked out in 
LORAN.) 



Number of lines per line group. 




Card 1, Format (40i2) 

Field 1, (Columns 1-2, 3-4, etc.), nu 

Number of lines in each line group to be treated individually. Must 
equal the number of line data cards read in for each group. When- 
ever a line is incorrectly assigned to a line group the program will stop 
and write out the message “LINE CENTER OUT OF GROUP FRE- 
QUENCY RANGE.” 


Data on line transitions . 


Cards 1,2. ..no. of lines, Format (i2,lx,i2,2el2.1,3el2.2,f5.1,f7.3) 
Fields 1-9, nd, nup, hvl, ff, gamp, xnol, gup, gl, el 
see descriptions below 


46 






Table DII. Index of Radiating Species and Level 


Index 

Species 

1-22 

N 

23-41 

0 

42-49 

C 

50-57 

H 

58-65 

Ar 

66-73 

N+ 

74-81 

0 + 

j 82 89 

C+ 


Field 1 (columns 1 2), nd index of the lower (absorbing) level of the line; identifies 
the radiating atomic and ionic species and level according to the scheme in 
table DII where the index increases (for each species) with increasing energy 
level 

Field 2 (columns 4-5), nup index of the upper (emitting) level of the line 

Field 3 (columns 6 17), hvl frequency of the center of the line, eV 

Field 4 (columns 18 29), ff f-number of the line when the lower level of the transi- 

tion is unlumped; when it is lumped, ff=(gee*f)/gee(sum over lumped), where 
gee is the statistical weight of the unlumped lower level of the transition 

Field 5 (columns 30-41), gamp HWHM per free electron due to Stark broadening, 
evaluated at 10 000 K (see ref. 3) 

Field 6 (columns 42 53), xnol number of lines in a given line group which have 

identical properties; it is used only when multiplet-averaged lines arc included; 
spectral properties are calculated for the averaged line; the line width is then 
multiplied by the number of original lines to obtain the effective total, only 
lines in the low-frequency region (roughly less than 5 eV) should be averaged in 
this manner; in cases where approximate results are sufficient, high-frequency 
lines may be multiplet-averaged as well 

Field 7 (columns 54 65), gup statistical weight of the upper level of the line (used in 
resonance broadening calculation only) 

Field 8 (columns 66 70), gl statistical weight of the lower level of the line 

Field 9 (columns 71 77), el energy level of the lower state of the line transition, eV 


Surface Absorptivity / Emissivity and Transmissivity Data 
Skip this input deck if kr(13) 1. 


Card Set 1, Free Format 

Fields 1 through nabtr, fabtr(j) 

Frequencies at which surface properties are given (eV) 
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Card Set 3, Free Format 


Fields 1 through nabtr, tmsw(j) 

Surface transmissivity at corresponding frequency 


NEQAIR Excitation and Radiation Inputs 

The following data are contained in file neq. in. They are read into the program in subroutines 
radin.f and molin.f. Note that all inputs are in cm - * for consistency with the NEQAIR code, 
but the inputs are immediately converted to electron volts within the code. 

Atomic Data 

This excitation database for atoms will not change unless the user has new data, new models, 
or new species. These data must be included for each atomic species in the radiation calculation. 

Energy level data. 



Card 1, Format (Ix,2a4,lx,gl0.3) 

Fields 1-2, (Columns 2-9), atomx 
Atomic species name. 

Field 3, (Columns 11-20), eiona 

Ionization potential of the corresponding ion (1/cm). 






Card set 3, Format (2i5,2gl0.3) 

Field 1, (Columns 1-5), me 

Index of the electronic level. 

Field 2, (Columns 6-10), npr 

Principal quantum number of the electronic level. 
Field 3, (Columns 11-20), elev 

Energy level of the electronic state (1/cm). 

Field 4, (Columns 21-30), stwt 

Degeneracy of the electronic state. 
> 

Transition data. 

( ^ 

Card set 1, Format (2i5,2gl 0.3) 

Field 1, (Columns 1-5), jl 
Index of lower level. 

Field 2, (Columns 6-10), ju 
Index of upper level. 

Field 3, (Columns 11-20), al 
Optical transition rate, 
s / 
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Card set 2, Format (il0,2gl0.3) 

Field 1, (Columns 1-10), jp 
Index of level. 

Field 2, (Columns 11-20), cz 

Constant multiplier in expression for continuum absorption coefficient. 
Field 3, (Columns 21-30), anz 

Exponent in expression for continuum absorption coefficient. 






Card set 3, Format (4(2i3,lpg8.1,0pf5.2)) 

The following fields repeat across and down this card set. 

Field 1, (Columns 1-3), jl 

Index of lower electronic level. 

Field 2, (Columns 4-6), kl 

Index of upper electronic level. 

Field 3, (Columns 7-14), al 

Constant multiplier in expression for excitation rate between two levels. 
Field 4, (Columns 15-19), el 

Exponent in expression for excitation rate between two levels. 

> 


Molecular Radiation Data 

The following inputs are used in the radiation calculation portion of the LORAN code and 
are the properties for radiating transitions. This database should be changed only if new data 
become available or new species are added. These data must be included for each molecular 
band in the radiation calculation. 
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Card 1, Format (Ix,2a3,3x,fl0.3,4i5) 

Field 1-2, (Columns 2-7), amol(l-2,ib) 

Name of band. 

Field 3, (Columns 11-20), aw(ib) 

Molecular weight of species, (g/mol). 

Field 4, (Columns 21-25), lmol(l,ib) 

Type index of molecular band. 

Field 5, (Columns 26-30), lmol(2,ib) 

Number of bands. 

Field 6, (Columns 31-35), lmol(3,ib) 

Number of energy levels. 

Field 7, (Columns 36-40), lmol(4,ib) 

Symmetry factor for heteronuclear (1) and homonuclear (2) molecules. 
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Card 2, Format (il0,7gl0.3) 

Field 1, (Columns 1-10), lmol(3,ib) 

As above 

Field 2, (Columns 11-20), sp(2,l,ib) 

Alternating line strength factor (not used). 

Field 3, (Columns 21-30), sp(3,l,ib) 

Dissociation energy of state (1/cm). 

Field 4, (Columns 31-40), sp(4,l,ib) 

Morse potential constant j3 (1/cm). 

Field 5, (Columns 41-50), sp(5,l,ib) 

Equilibrium internuclear separation distance (A). 

Field 6, (Columns 51-60), sp(6,l,ib) 

Dissociation energy including zero point vibrational energy (1/cm). 

Field 7, (Columns 61-70), sp(7,l,ib) 

A, a higher order energy level correction (not used) (1/cm). 

Field 8, (Columns 71-80), sp(8,l,ib) 

Angular momentum quantum number (not used). 


Card 3, Format (8gl0.3) 

Same as Card 2 for the upper state of the band (with 2nd index=2) except. 
Field 1, (Columns 1-10), sp(l,2,ib) 

Reduced mass of molecule (g/mol). 
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Card 4, Format (8gl0.3) 

Field 1, (Columns 1-10), sp(l,3,ib) 

Degeneracy of lower electronic state of band. 

Field 2, (Columns 11-20), sp(2,3,ib) 

Energy level of lower electronic state of band (1/cm). 

Field 3, (Columns 21-30), sp(3,3,ib) 

Spectroscopic constant ui e of lower state of band (1/cm). 

Field 4, (Columns 31-40), sp(4,3,ib) 

Spectroscopic constant u> e x e of lower state of band (1/cm). 

Field 5, (Columns 41-50), sp(5,3,ib) 

Spectroscopic constant u e t/ e of lower state of band (1/cm). 

Field 6, (Columns 51-60), sp(6,3,ib) 

Spectroscopic constant m e z e of lower state of band (1/cm). 

Field 7, (Columns 61-70), sp(7,3,ib) 

Spectroscopic constant b e of lower state of band (1/cm). 

Field 8, (Columns 71-80), sp(8,3,ib) 

Spectroscopic constant a e of lower state of band (1/cm). 


Card 5, Format (8gl0.3) 

Same as Card 4 for the upper state of the band (with 2nd index=4). 
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V 

Card 6, Format (2f5.1,2gl0.3,i5) 

Field 1, (Columns 1-5), vib (l,j,ib) 

Vibrational quantum number of upper state of vibrational transition. 

Field 2, (Columns 6-10), vib (2,j,ib) 

Vibrational quantum number of lower state of vibrational transition. 

Field 3, (Columns 11-20), vib (3,j,ib) 

Average value of the sum of the squares of the intermolecular distance 
r e for the vibrational transition. 

Field 4, (Columns 21-30), vib (4,j,ib) 

Franck-Condon factor for the vibrational transition. 

Field 5, (Columns 31-35), ka (j,ib) 

Maximum rotational quantum number for the vibrational transition. 

> 


Card 6 is repeated for each vibrational transition considered. 

M olecular Excitation Data 

The following inputs are used in the excitation calculation portion of the LORAN code and 
are the properties for all excited states required to obtain populations for those excited states 
which radiate. This database should be changed only if new data or species are available or if 
new excitation models require testing. These data must be included for each molecular species 
treated in the LORAN code. 
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Card 1, Format (Ix,2a3,3x,fl0.3,6i5) 

Fields 1-2, (Columns 2-7), amolec(l-2,is) 

Name of molecular species 

Field 3, (Columns 11-20), awb(is) 

Molecular weight of species. 

Field 4, (Columns 21-25), lmolec(l,is) 

Number of electron-impact dissociation cross-section data sets (maxi- 
mum 4 allowed). 

Field 5, (Columns 26-30), lmolec(2,is) 

Number of electron-impact, excitation cross-section datasets (maximum 
6 allowed). 

Field 6, (Columns 31-35), lmolec(3,is) 

Number of energy levels (maximum 8 allowed). 

Field 7, (Columns 36-40), lmolec(4,is) 

Number of first few electronic levels for which quasi-steady-state cal- 
culation is made (maximum 4 allowed). 

Field 8, (Columns 41-45), lmolec(5,is) 

Symmetry factor (l=hetero, 2=homonuclear). 
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Card set 2, Format (8gl0.3) 

Field 1, (Columns 1-10), spec2(l,jn,is) 

Electronic degeneracy of state jn. 

Field 2, (Columns 11-20), spec2(2,jn,is) 

Energy level of state jn (1/cm). 

Field 3, (Columns 21-30), spec2(3,jn,is) 

Spectroscopic constant u e of state jn (1/cm). 

Field 4, (Columns 31-40), spec2(4jn,is) 

Spectroscopic constant u e x e of state jn (1/cm). 
Field 5, (Columns 41-50), spec2(5,jn,is) 

Spectroscopic constant u e y e of state jn (1/cm). 
Field 6, (Columns 51-60), spec2(6 jn,is) 

Spectroscopic constant uj e z e of state jn (1/cm). 
Field 7, (Columns 61-70), spec2(7,jn,is) 

Spectroscopic constant b e of state jn (1/cm). 

Field 8, (Columns 71-80), spec2(8jn,is) 

Spectroscopic constant a e of state jn (1/cm). 

Field 9, (Columns 81-90), spec2(9jn,is) 

Equilibrium internuclear separation r e (Angstroms) 
Field 10, (Columns 91-100), spec2(10,jn,is) 

Dissociation energy £) 0 (1/cm). 

Field 11, (Columns 101-110), spec2(ll jn,is) 

Rotational spectroscopic constant D e (1/cm). 

Field 12, (Columns 111-120), spec2(12,jn,is) 

Second rotational spectroscopic constant p e (1/cm). 

S 

Card 2 is repeated for each energy level to be included for the molecule. 
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Card 3, Format (Ix,2a3,3x,7gl0.3) 

Field 1-2, (Columns 2-7), atom(l-2,l,is) 

First dissociated atom created from this molecular species. 
Field 3, (Columns 11-20), atw(l,l,is) 

Atomic weight of atom 
Field 4, (Columns 21-30), atw(2,l,is) 

Dissociation energy for the molecule (1/cm). 

Field 5, (Columns 31-40), atw(3,l,is) 

Degeneracy (statistical weight) of ground state of atom. 


Card 4, Format (Ix,2a3,3x,7gl0.3) 

The above information is repeated for the second dissociated atomic 
species. 


Card 5, Format (20a4) 

Field 1, (Columns 1-80), title 

Information concerning the source of the dissociation cross-sections 
which follow. 


Card 6, Format (8gl0.3) 

Field 1, (Columns 1-10), edis(it,is) 

Dissociation energy of it-th electronic state (1/cm). 
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Card 8, Format (8gl0.3) 

Field 1, (Columns 1-10, 11-20, etc.), sigmd 

Cross-sections for dissociation at energy ed (cm 2 ). 

Cards 5 8 are repeated for each set of dissociation cross sections. 

Card 9, Format (Ix,2a3,3x,3i5) 

Field 1, (Columns 2-4), spelll 

Lower state of molecular transition. 

Field 2, (Columns 5-7), spell2 

Upper state of molecular transition. 

Field 3-4, (Columns 11-20), ilst(l-2,it,is) 

Lower and upper electronic state index (ff=transition number in a 
molecule). 

Field 5, (Columns 21-25), ncard 

Number of cards of Franck-Condon factor data which follow. 
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Card set 10, Format (5(2i3,gl0.3)) 

The following fields repeat across and down this card set. 

Field 1, (Columns 1-3), jl 
Lower state index. 

Field 2, (Columns 4-6), j2 
Upper state index. 

Field 3, (Columns 7-16), fl 

Franck-Condon factor for the transition between states jl and j2. 


Card 11, Format (20a4) 

Field 1, (Columns 1-80), title 

Information concerning the source of the electron-impact excitation 
cross-sections which follow. 


Card 12, Format (8gl0.3) 

Field 1, (Columns 1-10, 11-20, etc.), ee 

Electron energy for electron-impact excitation cross-section (1/cm). 
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Card 14, Format (8gl0.3) 

Field 1, (Columns 1-10), arad(mu,ml,is) 

Radiative transition probability for the transition from level mu to ml 


Cards 9-14 are repeated for each transition considered. 
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Appendix E 
Tecplot File Format 


The Tecplot output produced by the LORAN code consists of both I- and IJ-ordered files. 
Each file starts with a variables declaration of the form shown in the first line of table EI. The 
second line contains information about the size of the file and the order of the data. In table EI, 
imax is the size of the I-ordered data set. If the optional jmax is included, it implies an IJ-ordered 
data set and gives a second dimension. The term order format is an optional variable indicating 
the data layout. When excluded, each column in the file corresponds to a variable. When set 
to block, each row of data corresponds to a variable (multiple rows will appear for each variable 
when the data set dimensions are large). Data for the zone then follow and are ordered according 
to orderformat and whether the file is I or IJ ordered. 

Additional zones of data can be included in the same file (for example, the variables along 
each normal grid line can be included as separate zones), as shown in table EI. 


Table EI. Tecplot File Format 


VARIABLES = vnamel, vname2, ... 

ZONE I = imax, [J = jmax,] [F = orderformat] 
Data for ordered zone 

ZONE I = imax, [J = jmax,] [F = orderformat] 
Data for ordered zone 

etc 
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